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LOW-STORAGE, EXPLICIT RUNGE-KUTTA SCHEMES FOR THE COMPRESSIBLE 

NAVIER-STOKES EQUATIONS 

CHRISTOPHER A. KENNEDY*, MARK H. CARPENTER* , AND R. MICHAEL LEWIS* 

Abstract. The derivation of low-storage, explicit Runge-Kutta (ERK) schemes lias been performed in the 
context, of integrating the compressible Navier-Stokes equations via direct numerical simulation. Optimiza- 
tion of ERK methods is done across the broad range of properties, such as stability and accuracy efficiency, 
linear and nonlinear stability, error control reliability, step change stability, and dissipation/dispersion accu- 
racy, subject to varying degrees of memory economization, following van der Houw r en and Wray, 1(5 ERK 
pairs are presented using from two to five registers of memory per equation, per grid point and having 
accuracies from third- to fifth-order. Methods have been assessed using the differential equation testing 
code DETEST, and with the ID wave equation. Two of the methods have been applied to the DNS of a 
compressible jet. as well as methane-air and hydrogen-air flames. Derived 3(2) and 4(3) pairs are competitive 
with existing full-storage methods. Although a substantial efficiency penalty accompanies use of two- and 
three-register, fifth-order methods, the best contemporary full-storage methods can be nearly matched while 
still saving two to three registers of memory. 

Key words, explicit Runge-Kutta, low-storage, numerical stability, error control 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. Direct numerical simulation (DNS) of the compressible Navier-Stokes equations is a 
means by which researchers may numerically probe the full range of scales in high-speed/fast-time-scale fluid 
behavior. Compressible DNS seeks to resolve all physically relevant, time and length scales associated with 
phenomena such as t urbulence, sound generation, and/or chemical reaction. Resolution of these phenomena 
is likely to require strict temporal error tolerances. The correspondingly accurate spatial discretizations 
involving possibly billions of grid points then fill the available memory of the computer. Hence, memory 
management of the time integrator is an important matter for DNS. The combination of high accuracy and 
low memory use potential for explicit Runge-Kutta (ERK) schemes makes them ideal for compressible DNS 
application. 

Efforts to reduce computer memory usage during numerical integration of ordinary differential equations 
(ODEs) have received sporadic attention in the past. [14, 27, 30, 56, 67, 90] For users confronted with severe 
computer storage constraints, established high-order methods such as the DOPRI5[21] may be prohibitively 
costly. Currently in the fluid dynamics community, users of ERK methods seeking to reduce memory usage 
have chosen to implement either a Williamson [87] scheme[25, 32, (54] or a van der Houwen[41. 42] (vdH) 
method. [4, 81] Williamson and vdH methods are both so called “2N” schemes, where N is the number of 
equations being integrated times the number of grid points. 

'Combustion Research Facility, Sandia National Laboratories, Livermore. California 94551-0969. E-mail: 
cak >§' boltzmann. ran .sandin.gov. 

t Aerodynamic and Acoustic Methods Branch, NASA Langley Research Center, Hampton, Virginia 23681-2199. E-mail: 
m.h.carpt n ter rc .nasa.gov. 

* Institute for Computer Applications in Science and Engineering (ICASE), NASA Langley Research Center. Hampton, 
Virginia 23681-2199. E-mail: buckarooKiicase.edu. This research was supported by the National Aeronautics and Space Admin- 
istration under NASA Contract No. NASl-97016 while the third author was in residence at ICASE, NASA Langley Research 
Center, Hampton, VA 23681-2199. 


1 



When solving the equation 

(II) ^ = F(1,V(t)) 

with an s-stage ERK method, a cavalier implementation requires the storage of the original ( /-vector, 
an intermediate ( —vector, and all s function evaluations. Williamson implicitly assumes that it is only 
necessary to be concerned with the memory requirement of a (/-vector, what is in effect a (/(/-vector, and 
that the memory requirement of F is inconsequential. Loosely, he implements the strategy over a single 
intermediate stage as 

(j dl ;U) = AjdUV-V + (A 

( r(j) = + BjdUU), 

where Aj and Bj are functions of the standard Butcher coefficients, F ij) and [ rU) are the j lh intermediate 
values of the function evaluation and the integration vector, and A1 is the step size. Note that (/, dU , and 
F must be stored. Unless work is done in a piecemeal fashion, three storage registers per variable will be 
required for the Williamson scheme. These methods have been referred to as 2N schemes. 

Wray [90], employing van der Houwen’s technique, places a set of conditions the scheme must satisfy 
that are more restrictive than W illiamson s. Van der Houwen and Wray devise a scheme where information 
alternates between the two available storage registers at each successive ERK stage. The procedure is loosely 
written over two intermediate stages as 

(Register 1) = X u) + (a j+lJ )AtF u) 

(Register 2) A' (/+1) = U u+1) + (bj - a j+lj )AtF ij} 

(>■ 3 ) 

(Register 2) i ■ r(j ' +2 > = A' (j+1) + a j+2 j +i AtF u+1) 

(Register 1) A' (j+2) = U u+ 2) + (6 j+1 - a j+ -, J+1 )AtF u+1) . 

By overwriting, the (', / . and A' vectors never fully coexist. The symbols a tJ and bj are the ordinary Butcher 
coefficients of the scheme. The .V vector may be thought of as a vehicle to bring information from previous 
stages into the current stage. To distinguish these methods from the Williamson class of 2N schemes, we 
will refer to them as 2R. schemes. 

f he primary difference in philosophy between the two methods is that in the vdff scheme, during the 
function evaluation, the previous solution vector is overwritten. Clearly there will be cases where this 
is not acceptable. Compressible DNS provides a situation where this method may be profitably utilized. 
This circumstance occurs because the ('—vector contains, principally, variables which are the products of 
variables needed to evaluate the flux terms. Consequently, the ('-vector must be decomposed into other 
variables, leaving the ('-vector itself disposable. Both Williamson (2N) and vdH (2R) schemes may be 
easily generalized to accommodate more than two storage registers (N or R). We make no claim that these 
two strategies are the only viable ones. We do suggest, however, that the vdH methodology is extremely 
aggressive in its conservation of computer memory usage. 

In the pursuit of computer memory use reduction, the first casualty is the retention of the (/'-vector at 
the beginning stage. Error control, in the more traditional sense, becomes impossible. A rejected step (such 
as violat ion of the error tolerance) cannot be restarted from ((<"> because with a 2N or 2R scheme, (/<”> is 
no longer available. Instead, alloting one additional register for an error estimate, one may monitor the error 
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occuring at each step and determine an appropriate next step. Including yet another register, ! could 
be retained so that action could be taken on a step which exceeds the predetermined error bound. This 
additional register approach, however, undermines the fundamental premise of this work and should not be 
used unless all other approaches fail. Schemes in this paper that, can be used in error monitoring/control 
mode are designated 2R+, etc., schemes because, if used in this mode, they require extra storage. Det ails 
of this implementation are contained in Appendix A. If overwriting is impossible, then the implementation 
must be modified and an extra storage register will be required. 

The goal of this paper is to derive broadly optimized, minimal-storage (2K+, 3R+, 4R+. and 5R+) 
ERK schemes based on only the vdH methodology and to explain how they are implemented. Choices 
are offered for storage reduced methods that address the needs of stability efficiency, accuracy efficiency, 
linear stability, nonlinear stability, dissipation and dispersion minimization, time-step/error control, and 
step-control stability. Invariably users will investigate physical phenomena that require different integrator 
properties, different error tolerances, and have different computer memory allocations. Hence, many good 
schemes are presented along with a rational basis by which to choose a scheme depending on one's needs. 
Based on the existing literature, the fluid dynamics community has been the largest customer for these 
low-storage schemes. For this reason (as well as personal research interests), optimization of ERK schemes is 
made with an eye towards the Navier-Stokes equations. Flows which are strongly viscous or chemically stiff 
may not be good candidates for these explicit methods. In addition, integrating the different ial-algebrah 
equations arising from the discretization of the incompressible or low Mach number equation set with an 
ODE-ERK method must be done with great caution so as to avoid drift-off and/or order reduction. 

A recurring criticism that accompanies use of high-order ERK schemes for discretized partial differential 
equations (PDEs) is the boundary “order reduction” phenomenon. [1, 13, 47, 49, 61, 66] Without proper 
care, order reduction occurring at the spatial boundaries can dominate the solution accuracy throughout 
the entire domain. The impact of the order reduction becomes more pronounced with increasing temporal 
accuracy. As such, the new schemes presented in this paper will be more susceptible to this problem than 
either Williamson’s or Wray’s third-order schemes. 

A second concern is that Runge-Kut.ta methods may seek out spurious fixed points of the differential 
equations being integrated. Methods exhibiting this behavior are called irregular. [2, 38, i8, 80, 83] All ERK 
methods greater than first-order accurate are irregular. W’e rely on error control and the fact that many 
equations are being integrated simultaneously to avoid spurious fixed points. 

Finally, we largely forsake aesthetic or “nice” coefficients (ones with simple, rational numbers) because 
the benefit from using a substantially more efficient integrator over hundreds of simulations, each taking 
tens or hundreds of hours, far outweighs the inconvenience of typing in twenty or so complicated coefficients 
correctly. Most solutions that are presented within this paper have been found numerically with established 
mathematical software. [18, 28, 29, 58, 88] Attempts were made t.o solve for schemes symbolically, but it was 
found that the assumption of various ay = 6j quickly made matters intractable because the equations of 
condition become algebraically nonlinear in the 6,'s. Scheme coefficients are given to at least 25 digits of 
accuracy. 

Some ERK background is necessary to facilitate a discussion on the optimization of accuracy efficiency, 
stability efficiency, error control reliability, step-control stability, linear stability, nonlinear stability, and 
dispersion and dissipation error within the context of storage reduction in later chapters. This will be done 
in sections 1 and 2. Two-register schemes will be reviewed in section 3 while three-, four-, and five-register 
schemes will be considered in sect ions 4, 5, and 6. Merits of the low-storage schemes are discussed in section 
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7, and comparisons are made with more traditional, full-storage ERK methods. In section 8, conclusions 
are drawn as to the utility of the various schemes. Appendices listing an implementation strategy and the 
relevant equations of constraint are also included. 

2. Background. We cannot hope to review the extensive body of Runge-Kutta literature germaine to 
integrating the Navier-Stokes equations. Therefore we tersely describe only those areas ofliterature that are 
crucial to the development of the new schemes being presented. For further details, appropriate references 
are provided. 

The compressible Navier-Stokes equations constitute a coupled set of partial differential equations that 
may be spatially discretized into a set of coupled ODEs with finite-difference techniques by the method of 
lines. We are concerned with the numerical solution of the initial value problem 

dU 

= F(tU(t)), U(a) = U 0 , /€M], 

where U = (p, pu, pc 0 , pYi) T is a function of the fluid density, p , velocity vector, u, total specific internal 
energy, c () , and species mass fraction, V' f . F contains the inviscid, viscous, reactive, and, possibly, body force 
terms of the compressible Navier-Stokes equations. 

Temporal discretization of the Navier-Stokes equations can be made with an s- stage ERK scheme, which 
may include an embedded error control scheme within the s-stage procedure. The implementation over a 
time Step A*, from time level t {n) to time level is accomplished as 

p{i) = inn = + * {i} = tw + aA /, 

(2.2) 

f ,( " +| ) = U {n) + At bjF u K fd»+ !) - {/(«) + At bjF ij \ 

wheie I —(()=;[ I) and = ( (/' n, -(-A^) are the solutions at time levels n and /i -f- 1 of order 

7 V T 1 and l ^ + ■ , the ( — vector associated with the embedded scheme, is of order p. The particular 
Butcher coefficients a,,, 6*. 6, , and c,- ol' the respective schemes are constrained by certain equations of 
condition, a short list ol which may be found in appendix B. In reading these conditions we remark that for 
a r/ ,h - order ERK, the k th equation of condition may be considered as [9, 35] 

(2.3) 7 ■<«> = -U<"> - 

a k q\ 

wlii< It defines < f > ^. , a scalar sum of Butcher coefficient products that, will appear throughout this paper. 
Both o and a vary with q and k. The conditions are identical to conditions with 6, replacing 6,. 

We assume that the standard row-sum condition applies: c, = a tJ . Extensive discussions of explicit 

Runge-Kutta methods may be found in the literature.[9, 21, 24, 35, 54, 69] Our style in this paper closely 
follows Dormand et al . [20] Schemes will be referred to as RKq(p)s[rR,nN+]X{ q dl5p , q djss }, where q is the 
order of the main scheme, p is the order of the embedded scheme, s is the number of stages, r/n is the 
number of registers used (vdH/ Williamson), + denotes that extra storage registers will be needed for error 
monitoring/control, X denotes either C (linear stability-error compromise), S (maximum linear stability), F 
(FSAL - first same as last), M (minimum truncation error), N (maximum nonlinear stability - contractive), or 
P (minimum phase error), and for “P methods, q d j sp and q d j ss are the respective dispersion and dissipation 
orders of accuracy. 
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2.1. Error and Error Control. Error in a </ th -order explicit Runge-Kutta scheme may be quantified 
in a general way by taking the L 2 and principal error norms, [ 02 , iSG] 


) _ || r (t/+l)| 


where r\ q) are the n q = { 1 , 1,2,4,9,20,48,115,286} error coefficients associated with order of accuracy 
q = { 1 , 2 , 3, 4 , 5. 6 , 7, 8 , 9}. For embedded schemes of accuracy p, additional definitions are useful, such as 

(2.5) 41 p + .) = || f (r + »||, - V (?J>' + ‘>) , 

N 3=1 

, (!1 

‘ 1 ~ * " l|f,,+,,|b ^ (r'i 

,,, ^ 


£■(»>+ 2 ) _ 


|f(?>+2) _ r (?>+2) 


D - Max{|«ij|, | 6 j|, 1 6,- 1 , |c,|}, 

p(,. + 2> _ ^ = = 

~~ ,A(r+D ||f(r +«)|| 2 


(r| ,,+2) ) 


One may also consider A& +1) , B { £+ 2) , d£ +2) , and E& +2) . All embedded schemes considered here are 
applied in local extrapolation mode; i.e., the solution is advanced with the higher - order formula. For 
a given order of accuracy, one strives to reduce A (<7+1) to as small number as possible. Both and 

C'(f+ 2 ) should be of order unity. The maximum magnitude of any of the Butcher coefficients, D< should 
be small, but may approach 20 in some high-quality pairs. [72] Shampine[<58] recommends £ (?,+2) < 1.5 and 
£ , 0 > + 2 ) < 0.5. Although these error measures are independent of the equations being integrated and hence 
onlv an approximate error metric, they will be used to select the "best scheme. \erner[ 8 o] points out that 
strictly relying on only # (/>+2) , and C (p+2) may not be adequate to distinguish among several good 

schemes. He also presents A& +1 \ A& +1 \ B { £ +2) , and C^ +J) . In another paper. Verner[84] recommends 
that should generally not vanish. Although not as frequently mentioned as the above parameters, the 

ratio of A^+ 2) /A {q+1) is sometimes controlled. For 5(4) pairs, Sharp and Smart[72] choose 5/2, Bogacki and 
Shampine[(j] limit it to 10, while Papakostas and Papageorgiou[(50] use 25. The risk of allowing .4 ^ /A 
to grow' too large is that the error controller may become less reliable at lax tolerances. Additionally, we 
require that all ^ 0 to avoid a defective embedded method, i.e., R : — 0. The stability domain of the 

embedded method is designed to be nearly as large as that of the main method to avoid instability in the 
embedded method at large step sizes. 

FSAL techniques, wdiere g s j — 67 , allow for the use of not only all function evaluations during an 
effectively s-stage computation, but also use + A After I (" + ') and + are evaluated, [A + ) is 
computed with s + 1 function evaluations. The principal motivation for doing this is that it allows more 
latitude in the design of the method and usually results in better schemes. The high stage numbers found 





in low-storage schemes make FSAL relatively less advantageous. Dense output via Runge-Kutta triples is 
forsaken here because there is little apparent interest within the DNS user community for such a feature. 
It may. however, find use if users seek global error estimates. [3] Pseudosymplectic or low-drift methods are 
also forsaken. 


2.2. Linear Stability. The stability function[36] for ERK methods is given bv 

(2. 10) «(c) = Det[<J, J -(«„ -e.ftjjr], 

(2.11) = 1 + *[ l) z + 4>' 2, .~ 2 + d>l 3) c 3 + + ^ 5) - 5 + . . . + 


with e _ {1,1. • ■ •. 1 ), S,j is the identity matrix, are the "tall trees,” and r contains information 
(eigenvalues) describing the equations being integrated. It is convenient in fluid dynamics to consider linear 
stability in the context of the prototypical, one-dimensional, convection-diffusion equation 


(2.12) 


dU 

Hi 



d 2 

+ < *'ae 



where a is a convection or sound velocity and a v is a diffusivity[48] (mass, momentum, or energy). Other 
st udies of st ability of ERK methods applied to the compressible Navier-Stokes equat ions have been conducted 
by Sowa[75] and Muller. [59] If the spatial derivatives are considered as high-order, centered, finite-difference 
operators then the Fourier image of the convection-diffusion equation becomes 

(2.13) r = -A* + A, 

(2 14 ) xjr — 'pa sin(£) + 2 b sin(2g) + 2 c sin(3£) + • • •] 

[1 -f 2a cos(£) + 2/3 cos(2£) -f ■ • •] 


In this expression A — and A„ = are the inviscid and viscous CFL numbers, Ax is the local spatial 

grid spacing. At is the magnit ude of one time step, 0 < £ < 7r is the spatial wavenumber, 4* is the Fourier 
image of the first, derivative operator, and {a, b, c, a, /?} are coefficients of the derivative operator used in 
evaluating the convection-diffusion equation. As the “compact/’ sixth-order derivative operator is popular in 
compressible DNS, these last coefficients will be set to {7/9, 1/36, 0, 1/3, 0 }.[48] 

A stable method has |/?(r)| < 1 at a particular value of r or for ail wavenumbers at the pair (A,A r ). 
This requirement is necessary but may not be sufficient, [33, 50, 65] Unlike many contemporary ERK pairs, 
imaginary axis stability is a high priority to the methods designed in this paper. The derived linear sta- 
bility domains, in terms of (A, A*,), are a strong function of which derivative operator is chosen. Reducing 
truncation error of the spatial derivative operator reduces the extent of the linear stability regime. Use of 
the corresponding second derivative operator rather than repeated use of the first derivative operators for 
the viscous terms reduces the maximum viscous CFL number. Nevertheless, determining linear stability 
as previously described gives results representative of a broad class of numerical methods used for DNS of 
compressible flows. 


2.3. Step Control Stability. We consider two step-size control strategies: [34, 36, 37] 


(2.15) 

(A f) ( " +1, = K(A<) 

(2.16) 

(A<) (w+1) = k(A ty 


(n) 


_ - r 

P" +i) ikj ■ 


where e is some chosen integration error tolerance, k & 0.9, and <W 1+ i) pr(» + i) _ £/0* + 1) an( j 

most common method, Eq. (2.15) is an example of an integral feedback (1-) controller. The second, more 
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sophisticated, Eq. (2.16) adds a proportional feedback component and is called a Pi-controller. Following 
Hairer and Wanner, [36] we define 


(2.17) 

<t 

R(:)= E(z 

)= E 



/. = 3f 

[«'(--)-] 
. RU) . 

. V 

= 


i=l 

>=/’+ 1 







as well 

as the matrices 












( 1 

u 

0 

0 

\ 


( i 

« \ 


—a 

(1-t 

r> u ) 3 

3r 


(2.18) 

6 “ ( (i 

- an) ) ’ 

c = 

1 

0 

0 

0 







1 

0 

0 

/ 


In the case of the C-matrix corresponding to Eq. (2.15), a = (p+1) For C corresponding to Eq. (2.16), we 
set. q = 0.7 fp and 3 = 0.4/p where p is the order of the embedded method. If at the regions where |/?(~)| = 1 
the spectral radius of C or C is less than unity, t hen the step-size control mechanism is said to be SC-stable. 
As DNS runs are often made near the linear stability limits of the integrator, step-size change oscillations 
may result and give rise to a rapid accumulation of global error during oversteps. Still more involved would 
be a PID-controller, which we do not use. The PID-controller is obtained by multiplying the RHS of Eq. 
(2.16) by {c/||rf (n ~ 1) ||ocP and creating a 6 x 6 matrix C which has the elements C’.j = Qj,iJ = 1,2, 3, 4, 
(% 5 = — (\ c = c, and CV ,3 = Cf ,4 = 1. with all remaining elements being zero. 

2.4. Dispersion and Dissipation Error. Dispersion and dissipation of ER.K methods[10, 44, 4o] may 
be considered by taking the derivative of V = c twl with respect to time, dU/dt = iulU , where u> is a temporal 
wavenumber (frequency). The stability function for this ODE has the argument w where v — ui(At). An 
ERK method where R and I, respectively, are the real and imaginary parts of R{iu) is said to be dispersive 
of order q<j,sp and dissipative of order q<a ss if 


(2.19) <£(*/) = v — arg (R(iv)) = v— arctan 


n 


2j + 1 




Qdisp 4" 1 ] 


j=» 


(2.20) o(i/) = 1 - \R{iu)\ = 1 - y/ft- + T 1 = = °( t/<,d ' ss + 1 )- 

i=o 


Hence R(w) = 7v . + il = (1 — a(t/)) f •(•'-♦l 1 ')). Some authors refer to the phase-lag order of the method, 
which, in our notation, would be (gdisp.tfdiss - 1). Control of both spatial and temporal dissipation and 
dispersion in acoustics applications has employed ERK methods satisfying only quadrature[46, 92] and 
subquadrature[93] order conditions. Applied to nonlinear problems like the Navier-Stokes equations, these 
methods will generally not exceed second- and third-order accuracy, respectively. 

2.5. Nonlinear Stability. Nonlinear stability of Runge-Kut.ta methods[53] focuses on the discrete 
analog to dissipativity of F(t,U(t)) in some given norm, 


( 2 . 21 ) 


|| U(t + At) - U(t + A/) 1 1 < || U(t) - C(0l|. 


where V is a perturbed approximation to V and F(t,V(t)) belongs to any one of the four function classes: 
linear (£) or nonlinear [T] and dissipative in an inner product or maximum norm. For ERK methods, [15, 52, 
53, 55, 91] the dissipativity criterion is replaced with the so-called circle-condition, and maximum step-sizes 
are related to a contractivity threshold: the largest possible step-size that ensures 1 1 / f ” + 1 _ F ! ” + 11 1 1 < 



IK /(n) - T'” , ||. A radius of conditional or circular contractivitv for the four function classes may be denoted 

rr 3 , r £j where r^< r£<x < r c ,, and r^< r Ti < r Cj . We will call a method (conditionally) 
contractive if at least, vjr 2 > 0. 

K raaije vanger [52] (C^,) and Dahlquist and Jeltsch[16] (? 2 ) have shown that no ERK method has r T > 0 
and is greater than fourth-order accurate. Maximum norm contractivity is closely related to positivity. [39] 
Positivity is particularly appealing because it ensures that physical quantities such as temperature and 
species concentrations remain forever nonnegative. The radius of positivity is the same as the radius of 
maximum norm contractivity. 


To determine we first define the mat riees[16] M (j = 6,-a.j + b jaji - b { bj, Bij = diag {6] 6, } , 

and Mu — B lJ ' . If bj > 0. then ry_, = —jjbr where is the smallest eigenvalue of the 

matrix Mu. A nonvanishing value of ^requires that b, > 0, a {j > 0, and the Runge-Kutta K-function, 
K(Z) = det[/ - (A - eb r )Z], is absolutely monotonic on [-r^, 0] where Z = {r,, to, • • • , r*}. 1\(Z) is said 
to be absolutely monotonic at a point £ if[5'2] 


(2.22) 


d‘ l +ij+ +it K{£,£, ■■■,£) 
d;Vdzi 3 ---d& 


> 0 . 


In the case of ERKs, each ij may be equal to either 0 or 1. The largest magnitude of £ on the negative 
real - axis for which these 2* inequalities hold is denoted -r Toe . Alternatively for ERKs, one may enforce 
nonnegativity of 


(2.23) R(£), A(£) = A(1 - £A)~\ B(£) = b T (I -*A)-\ £(£) = (I - ^ A) _1 e, 

at s = where e = {1, 1, ■, 1}, b = 6,, and A = a ,, . Assuming that b, > 0, a,, > 0, these present 

1, (s - 2)(.s - 1 )/2, (s - 1), and (s - 1) inequalities, respectively, or s(s + ])/2 in total. It should be noted 
that the region of circular contractivity is a circle located at z = -r with radius r r ^_ or , whichever is 
appropriate. 'I his implies the property is likely most useful for parabolic rather than hyperbolic equations. 
For comparison purposes we follow Dahlquist and Jeltsch and write r £ , = , the corresponding radius of 

the linear problem in an inner product norm, [57] i.e., the largest circle centered on the negative real - axis, 
fully contained in the left half-plane, that fits within the region where |/f(;)| = 1. The linear analog of 
r T * c > s »•<.% • The stability function, R(Z), is said to be absolutely monotonic at a point £ if[53] d’R(£)/d:‘ 
> 0, i = 0, 1.2, • • -. 8 . The largest, magnitude of £ on the negative real - axis for which all of these *• + I 
inequalities hold is denoted —r £ae . Kraaijevanger[51] gives the maximum achievable r £oo per stage for an 
m-stage method with order his optimal scaled threshold factors. We do not consider the internal stability 
of ERKs. [42, 43] Nonlinear instability caused by spurious triad wave interactions[79] in the spatial domain 
is outside the scope of this paper, but is probably best dealt with by using high-order filtering. [48] 

2.6. Efficiency. Efficiency of a given s-stage ERK scheme may be considered from two decidedly 
different perspectives. One philosophy assumes that temporal integration error is acceptable and seeks to 
time step as briskly as possible. Simulations running on expensive supercomputers for hundreds of hours are 
under great pressure to be integrated as quickly as possible. Alternatively, integration may be conducted at 
some chosen maximum acceptable error. Virtually all DNS efforts that, these authors are aware of implicitly 
subscribe to the former philosophy, due to computer resource limitations. 

Stability-limited time stepping is the more primitive approach and only seeks to compare the relative 
efficiency of two schemes by using[69] 


(2.24) 


j (st.ab) Ai / ^y 1 = (Af>, 

' S 1 V *‘2 / A 1 
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where A is understood to be either the inviscid or viscous CFL number and scheme 1 is most efficient for 
;/ (stab) > j q-jjjg term compares the distance integrated per unit of work (evaluations of F(t. U(1)) — number 
of stages) with no regard for the accuracy of integrat ion and may be used to compare met hods with arbitrary 
orders of accuracy and numbers of stages. For viscously or reactivelv dominated problems this term could 
be amended by replacing the CFL numbers with the respective ^ T 1 or °f each scheme. 

Relative efficiencies of two (/"'-order schemes based on an error limited time stepping procedure might 
best be measured by [40, 09] 


(2.25) 


(acc) 


' Ai, q+l) ' 


*i \A 


(9 + 1 ) 
1 


(AQt « 3 

(A/)j ’ 


where scheme 1 is most efficient for i^ cr) > 1. Slightly different from (/ st:lb) , r/ (acc) compares the distance 
integrated per unit of work at fixed integration error. (A/)’. We will consider the number of stages in a 
FSAL method as the effective number of stages for efficiency purposes. Obviously, for sufficiently large error 
tolerances, large time steps might exceed the linear stability bounds. In comparing schemes with different 
orders of accuracy we cannot use this last metric and simply paraphrase Shampine[69] . For sufficiently small 
error tolerances the higher order method is more efficient, but this argument does not imply that, the lower 
order method is more efficient for large error tolerances. Prince and Dormand[62] note that only on the 
linear problem are lower order formulae sometimes preferable to higher order formulae. Sharp[71] finds that 
the higher-order methods are generally more accuracy efficient on nonstifT equations. 

The choice of which efficiency measure should be used depends most strongly on what level of error is 
acceptable to the user. This, in turn, depends on what physical phenomena are being sought through the 
calculation. If integration at the linear stability limits produces sufficiently small error then, efficiency is best 
considered by using >/ (stab) ; otherwise i/ (acc) seems more appropriate. Spatial accuracy, another important 
matter that we do not consider here, must also be addressed. Strict temporal error tolerances make little 
sense without correspondingly strict spatial error tolerances. A future study on the spatial and temporal 
accuracy/resolution requirements associated with particular physical phenomena would be of tremendous 
utility to compressible DNS practitioners. 


2.7. Simplifying Assumptions. Finally, in the course of designing several of the schemes in this 
paper, resort, ing to Rut,cher[8, 3(3] simplifying assumptions will be useful. On occasion, assumptions 


(2-26) 

C(n) : 

H a *7 c l 1 = 7' *‘ =1 s > q ~ v 11 

j = l 1 

(2-27) 

D(Q : 

II 

I 

<ri > 
ii 

e 

-Wii 


will be invoked. 


o 



3. Two-Register Schemes. An s-stage ERI\ method placed in two-register vdH format, (see van der 
Houwen[42], equation 2.2.4') takes on the Butcher array form 



and allows (2s — 1) degrees of freedom (DOF) to satisfy all constraints. In general, for an r-register method, 
there will be r * s - r * (r - l)/2 DOF available. Conversely, (s - r + 1) * (s - »■)/ 2 DOF are sacrificed for 
low storage. Setting r = s, the basic ERK method is retrieved with s * (s + l)/2 DOF. 

3.1. Two-Stage, Second-Order: RK2()2[2R]. All two-stage, second-order ERKs may be used in 
2R format. Minimum A^ 3) = 1/6 occurs at c 2 = a 2] = 2/3, 6j = 1/4, b-, = 3/4 with (r^, = 

(0./91,0.500). The maximally contractive second-order method is Huens method;[54] c-> = «•>] = 1, 
— ^ ■ ^2 = l/'2 where Vr 2 = = 1 and T' ! * = \/2/6 - These are of only academic interest to the 

compressible DNS community because when implemented with centered, finite-difference methods on the 
convection-diffusion equation, the methods are unconditionally unstable in the inviscid limit. 

Tin ee- Stage, Third-Order: RK3(2)3[2R-j-]. The general solution to the two- register, three- 
stage, third-order vdH scheme has a one-parameter family of solutions along with two specific cases. These 
are readily derived by following Lambert ’s[54] three cases. For the cases when c 2 ^ (),§, or c 3 and when 
^ 0, the one parameter, C3, family of solutions is obtained for 


4 — Tea + 6c§ ± y/c$ ( 17 — 60c 3 + 84cf( — 48r:|) 
6(1 -2c 3 + 2ci) 


Provided c 2 is not complex. From this, Wray[90] suggests c 3 = 2/3, yielding c 2 -- 8/15. Minimum principal 
error norm for the RK3(2)3[2R+]M is found by solving 0A i4) /dc 3 = 0 for c 3 by using the minus solution 
above. The result, is that the minima occurs at c 3 ss ^0777641 where A {4) = 0.04412 and (r^„, vyr^) 
(0.-)21. 0.150). Maximum vyr ^ occurs with RK 3 ( 2 ) 3 [2 R -)-] N at c 3 ss fj'7703 by using the plus solution where 
,4 (li = 0.05094 and ( vjr , . = (1.127,0.838). Asking for contractivity of the embedded methods had the 

unfortunate consequence of increasing £ ,(4) above optimal. The minimum principal error norm achievable for 
any explicit RK3()3 is ,4 (4) = 0.011809. Maximum radius of contractivity for the general RK3()3 is = 1 
or r_f 2 sa 1.215. All RI\3()3 methods have »y ^ = 1 and = 1.256. 

In the two specific cases where r 2 = c 3 = 2/3 and b 3 = 3/5. then A w = 0.04630 and (rjr,, r ? x ) 
= (1.0, 0.6), and where c 2 = 2/3, r 3 = 0, 63 = (1 ^ \ZV7)/S, both solutions have ,-l !4 ' = 0.1326 and are 
noncontract ive. The former confluent solution admits only a defective embedded method. Stability limits of 
all three-stage, third-order, ERK schemes are (A, A„) = (0.87,0.63) when integrating the convection-diffusion 
equation discretized with a sixth-order, tridiagonal, first-derivative operator. 


3.3. Four-Stage, Third-Order: RK3(2)4[2R-f]. Using two-registers over four stages allots degrees 
of freedom. Enforcing third-order accuracy, r<*> = 0, £ = 1,2,3, leaves three remaining DOF. For accuracy 
efficiency, RI\3(2)4[2R+]C minimizes .4 (4 ) subject to = 1/24 in order to maximize linear stability 
and dissipation order. The resulting scheme is 6% more accuracy efficient, than RK3(2)3[2R+]M, and has 
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(A. A, ) = (1.42,0.70); this scheme is listed in Tables 1 and 4 and shown in Figure 3.1. For nonlinear 
stability, RK3(2)4[2R+]CN seeks maximal while achieving sixth-order dispersion error, <*> 6 = 0. by setting 
<j)^l = Of the maximum possible rj: cc = 2 for any RK3()4[52] with A 1 ' 1 1 = 0.03608, RK3(2)4[2RT]C N 

achieves r^= 1.007 while keeping A™ = 0.02870. As with the RK3(2)3[2R+]N. a contractive embedded 
method drove £- 4) to slightly greater than 1. If < f!/ ) = 1/48, giving (A, A,.) = (1.08, 1.30), 2 is ]>ossible 

for RK3Q4 methods. 


3.4. Five-Stage. Fourth-Order: RK4(3)5[2R+]. Adding a fifth stage to a 2R-vdH scheme provides 
nine degrees of freedom. Fourth-order accuracy may now be considered. Eight order-of-accuracy constraints, 
T (k) — o j. — i 2, 3, 4, leave one DOF to optimize linear stability while maintaining acceptable accuracy via 
variation of 4>i/\ Tables 1 and 3 and Figure 3.1 give the $!, * = 1/206 solution, RK4(3)5[2R.+]C, with 
.4(5) = 0.005121 and (A, A, ) = (1.67,1.21). Mated to this is an embedded scheme with 4 >( 1 ‘ l 1 = 1/28. We 
were unable to find any contractive methods for the RK4(3)5[2R+] or phase-lag methods having reasonable 
principal error norms. Setting < I>i l 5 * = 1/240 gives (A, A,,) = (1.62, 1.48) and the largest rf^for the RK4()-> 
methods is 2. 

3.5. Six-Stage. Fourth-Order: RK4(3)6[2R+]. As addit ional stages can somet imes make for more 

efficient methods, [72] one may consider an RK4(3)6[2R.+] scheme with three residual DO! aft.ei satisfying 
r (fc) —oi;— i 2, 3. 4. Searching for solutions uncovered RK4(3)6[2R+]C with ,4 (o ) = 0.002148 and (A, A,, ) = 
(1.97,1.18). Unfortunately both r/ acc) and r/ (stab) are less than those of RK4(3)5[2R+]C. No attempt was 
made to find a contractive solution. At $J, 5) fa 1/159 and » 1/2529 where (A, A, ,) = ( 1 .85, 1 .62), r £ao may- 
reach fa 2.651. For increased phase-lag accuracy one may set <j>r, = <J>7 = 0 to find (A, A,.) = (0.29,0.96) or 
may set 0 5 = q 6 = 0 to find (A, A,,) = (0.35,0.89). Minimizing dissipation error with o 6 = a 8 = 0 results 
in $<, 5) = 1/128, = 1/1152, and (A,A„) = (1.94, 1.09). With A w = 0.002509, RK4(3)6[2R+]P{4,9} is 

such a scheme. Both RK4(3)6[2R+]C and R.K4(3)6[2R+]P{4,9} use an embedded method with t!, 1 ’ = 1/26 
and $!, r,) = 1/150. 

3.6. Nine-Stage. Fifth-Order: RK5(4)9[2R+]. A fifth-order, 2R-vdH scheme may be obtained in 

nine stages by solving the 17 unsimplified equations of condition, r (fc) = (),!• = 1,2, • • •,5, for the 17 free 
Butcher coefficients. Solution properties cannot be optimized. Over 800 distinct real roots to this system 
of equations have been found. The most accurate root found. RK5(4)9[2R+]M, has -4 (l,) = 0.0006172, but 
(A, A„ ) = (0.21 , 1.03). The most stable method, RK5(4)9[2R+]S, has .4 ' * 1 = 0.001014, ( A, A, ) — ( 1 .78, 1 .o9), 
and (A, A„ ) = (1.60.1.61). A compromise solution, RK5(4)9[2R+]C, was found with .4 (6) = 0.0008209. 
(A, A, ) = (1.05. 1.29), and (A, A, ) = (1.63, 1.15). An embedded method was designed for the three schemes 
by satisfying all eight fourth-order constraints plus setting $j, 5) = 1/135. These methods are presented in 
Tables 1 and 6. Stability diagrams are provided in Figure 3.1. The largest linear positivity radius for these 
RK5(4)9 methods appears to be fa 4.095 occuring at } 1/779, $!,* fa 1/7444. $ ( *!j fa 1/121935, 

fa 1/4494000, where (A, A„ ) = (0.34,2.24). 

4. Three-Register Schemes. Applications having slightly less stringent memory constraints may add 
an additional storage register per ODE. Extending the vdll methodology to three-registers, an s-stage scheme 


li 



takes the Butcher array form 


o 


c> a 21 
O a 32 

6] a^i CI43 

f >2 U 53 (1 54 


where there are now (3 
Alternatively, (s — 2) * ( 


£><-3 a $,$—2 a 5 , 5—1 

^2 **' t>s— 3 6 s _2 6 s _ 1 6 ? 

s — 3) independent coefficients that may be used 
s - 3)/2 coefficients are lost to low storage. 


to satisfy particular 


conditions. 


4.1. Three-Stage, Third-Order: RK3()3[3R]. Any three-stage, third-order ERK method maybe 

implemented in 3R format. As such, one may seek the method having the smallest principal error norm. 
By setting dA^/dc? = 0 and dA^ /dc 3 = 0 from the two-parameter family of solutions it is found that 
r? ~ M 79 5 282 5 ' C 3 ^ 326524 — 0.041809, and v t^) — (0.894, 0.0). The relation of c 2 and C3 to the 

various other Butcher coefficients may be found in the literature. [9, 21, 35, 54] Stability limits are identical 
to the RK3(2)3[2R+] methods. Maximal contractivitv, r? oo = r?, — 1, is found in Fehlberg s[26, 52, 74] 
method with c 2 = 1, C3 = 1/2, and = 0.07217, while for Cooper s scheme[15] in an inner product norm 
where (r^ 3 , r^J = (1.215,0.691), r 2 * 270/251, c 3 % 166/305, and A (4) = 0.07221. 

4.2. Four-Stage, Third-Order: RK3(2)4[3R+]. Kraaijevanger[52] has shown that optimizing the 

radius of maximum norm contractivity for general third-order ERKs allows one to obtain (s - 2) 

for a = 3,4. tor s > 5, < (s — >/s). A family of third-order schemes given b} r four unique Butcher 

coefficients, 6, = (s - 2)/(.s(s - 1)). * = 1, 2, • • ■ , (s - 1), b 9 = 2/s, a,, = l/(s-2), i = 2,3, • • 1) > j, 

<hj = l/(2(«- 1)), i - 1 , 2, • • , (s - 1), which for s = 3, 4 constitute the maximally contractive 

methods. From these relations it is seen that c,-, = 1 and c, = 1/2. For s - 5, 6, 7, and 8 one finds for 

this family that (ryr 3 , vyr ^ ) is given by (2.449,2.202), (2.828,2.347), (3.162,2.460), and (3.464,2.553). For 

reduced storage we set = (s - 2)/(#(* - 1)) = « s , = l/(2(s - 1)) to find ,s = 4. The resulting method. 
RK3(2)4[3R+]N. is essentially given by Kraaijevanger with r T , = >'c^ = ' r . = 2 and .-1 (4) = 0.03608. 

A good embedded method for this scheme is bj = {8, 9, 8, 60J/85. 

4.3. Four-Stage, Fourth-Order: RK4()4[3R]. From the general solution to the four-stage, fourth- 
order ERK scheme, [9, 21, 35] it is found that there is a one-parameter family of 3R solutions and three 
specific solut ions. The one-parameter family of solutions is given by 

, 4 n _ (20 - 50c 2 4- 36<-g) ± ^ (-20 + 50c 3 - 36ciQ 2 - 4(9 - 26c 2 -I- 16e-g)(16 - 36r-> + 36c;l 
3 2(lfi-36r 2 + 36c2) 

where 0,e 2 ,c 3 , 1 are all distinct, c 2 # 1/2, 3 - 4(c 2 - c 3 ) + 6c 2 c 3 ^ 0, and r 3 is not complex. The principal 
error norm is minimized by setting dA {5) /dc 2 = 0 where RK4()4[3R]M is found for the plus solution with 
c - ^ lsuncmii ' rjr -j = 0.718, and ,4* 5 ' = 0.01263. Maximal r? 3 = 0.882 occurs with the plus solution 
R.K4()4[3R]N at r 2 as (with .4< 5) = 0.01319) for 3R methods and r? 3 = 1.144 for Cooper’s [15] R.K4()4 
method. These values compare with r 7i = 1 and ,4 (5) = 0.01450 for the “classical Runge-Kutta” (see 
Butcher[9], §313) and .4 (o) = 0.011977 for the absolute minimum principal error norm for any four-stage, 
fourth-order, ERK scheme. Kraaijevanger[52] has shown that there exists no L 0 c contractive RK4()4 method. 
Adding a third-order embedded scheme to this method is impossible unless FSAL techniques are used but 
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we do not pursue this matter. Instead of a FSAL pair, complete use of the fifth stage generally makes for 
more efficient schemes. One potential exception is inviscid stability efficiency (A/* « 0.355). The three other 
specific cases are found using Butcher's eases 3, 4, and 5 in §312. [9] In case 3, by = —2/15 and .4^0 = 0.03416. 
in case 4, b 4 = 3/10 and ,4 (5> = 0.02330, while in case 5, c 2 = 3/7 and A {:,) = 0.01282. None of these last, 
three schemes is contractive. Linear stability limits on the convection-diffusion equation for the RK4()4[3R] 
are (A, A„) = (1.42,0.70), r £eo = 1, and r £ , = 1.393. 

4.4. Five-Stage, Fourth-Order: RK4(3)5[3R+]. Three additional degrees of freedom afforded 
by adding a third register to the RK4(3)5[2R+] method may be put, to good use. Optimizing accuracy, 
RK4(3)5[3R+]M has T (5) = 0.001884 and (A.A„) = (0.22,0.81) where tv$ = 0. A similar method of 
Prince has r]^ 89 = 0. A balance between linear stability and low error is found in RI\4(3)5[3R+]C with 
.4(5) = 0.003859 and (A, A,.) = (1.67,1.17). It. should be noted that for the RK4(3)5[3R+]C, selection of 
= 1/200 forces A< 5) > 0.003333, as can be seen from .4^. Contractivit.y appears to be maximized with 
RK4(3)5[3R+]N having (r^ 3 , vjr^) =(0.995,0.477), 4 (5) = 0.004587, and (A. A,.) = (1.67, 1.20). Although 
not nearly as contractive as Kraaijevanger’s RK4()5 scheme, it has 7% better »/ acc *. Each of these three 
schemes is presented in Tables 2 and 5. Stability plots are given in Figure 4.1. Two highly accurate 
RK4(3)5[3R+]P{4, 7} schemes where 4>|, r,) = 1/144 and (A, A,,) = ( 1.74, 0.89) were found with .4 (5) = 0.002658 
and A (5) = 0.002857, but. neither would accept an embedded method with a reasonably large linear stability 
region . 


4.5. Six-Stage, Fifth-Order: RK5()6[3R]. With only 15 degrees of freedom simplifying assumption 
D(l) may be invoked to reduce the number of condition equations from 17 to 15. By doing so, a fourth- 
order embedded scheme is no longer possible. At least 13 schemes like this exist.; the most, accurate found. 
RI\5()6[3R]M, has .4 (fi) = 0.003678 with (A, A,,) = (0.20,0.72). 


4.6. Seven-Stage, Fifth-Order: RK5(4)7[3R+]. To get. a 5(4) pair, a seventh stage is added and 
only simplifying assumption 0(2) is utilized. This results in 18 equations in 18 unknowns foi the main 
scheme and 7 equations in 7 unknowns for the embedded method. Of the 7 schemes found, the best one is 
R.K5(4)7[3R.+]M with ,4< 6 > = 0.002213, (A. A,) = (0.28,0.92), (A, A,) = (0.95.0.59), and B < 6 > < 1.0. Adding 

lead to a method with as much as 38% better as will 


an extra stage to this method, however, ran 
be shown in section 4.7. Maximum 2.(554 occurs at. 

(0.28, 1.60). 


(G) 

20 


1/955. 


]/l 7715)5 with (A. A r ) = 


4.7. Eight-Stage, Fifth-Order -RK5(4)8[3R+]. An eight-stage, three-register vdH scheme has 21 
degrees of freedom. Seeking a 5(4) pair, Butcher simplifying assumption (7(2) is applied. The resulting 
system of equations necessary to satisfy all order conditions is 


(4.2) 


Jfc) 


J4) 


: 0 , 

( 5) 

r 4.5,8 


k = 1, 2, • • • , 5, dijCj = cj /2, ? = 3, 4, ■ • * , 8, 

= &2 — biClj? = 5Z,-_3 hi C{(li2 — E«j=3 6 »* a «J°j2 “ 


for the main scheme and 

5 

(4.3) r 1 (/r) = 0, k = 1,2, 3, 4, 6 2 = r] 4) = M 12 = 0, 

for the embedded scheme. Optimization may now be done with two remaining DOF in the main method 
and one in the embedded method. A numerical search found two low-error solution families (among 25 or 
so), the first with more desirable stability properties and the second having lower ^ • RKo(4)8[3R-b]( . 
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RI\5(4)8[3R+]P{8,7}, and RK5(4)8[3R+]M are given in Tables 2 and 6. The first two come from the more 
stable family. RK5(4)8[3R+]C has .4< 6 > = 0.0008306 with (A, A„) = (1.30, 1.52) while RK5(4)8[3R+]P{8 ,7} 
has .10') = 0.0007923 with (A, A,) = (1.01, 1.20). RI<5(4)8[3R+]M achieves .4< c > = 0.0003240, but (A,A„) = 
(0.32, 1.00). The final degree of freedom for the embedded methods is used to set T,*, 5 * = 1/130, 1/135, 1 /122.5 

in the RK5(4)8[3R+]C. P, M schemes, respectively. Stability plots for the three schemes are shown in Figure 
4.1. 


Enhanced dispersion/dissipation order is enforced with 


*20 = (1 + 22684> l ”' 5 )/756 <J>^ = (1 + 22680<F ( , 8 1 ) 5 )/7560 (0 7 = <p 9 = 0), 

= ( 1 + 57604> , 1 8) 5 )/5760 (q 6 = o 8 = 0), 

1/5040 (o 6 = 0 7 = O). 


si**) \ /TC 


h< 7 > - 


(4.4) 


= 1 /720 


* { m ] = 1/720 


4> 


48 


With the RK5(4)8[3R+]M solution family, the methods RK5(4)8[3R+]PM{10,5), RK5(4)8[3R+]PM{8,7}, 
and RK5(4)8[3R+]PM{6,9} may be found having A <6) = 0.0005049, 0.0005946, and 0.0005525, and (A, A„) = 
(0.35,1.14), (0.88,0.98). and (0.53,1.04), respectively. Each may be fitted with a high-quality embedded 
method by setting t!, 5 ' = 1/130. In each of these methods, D < 2, /#'•'’) , C fi ) < J.5 and T 1 * 6 ) < 0.5. The 
largest possible r^for RK5(4)8 schemes is found at < t >[^ rs 1/834, 1/9862, rs 1/266413 where 

3.368 and (A, A,,) = (0.30, 1.89). 

5. Four-Register Schemes. Further relaxing the memory constraints, the 4R-vdH scheme structure 
appears as 

o 

I 

a 2l 

Cl 3 | 

a 4 1 


a 3 2 

(14 2 


«43 


«52 «53 


3 


«54 


a 64 




«65 


^1 ^2 b s 

and has (4s - G) DOF. Storage reduction has consumed (s - 3) * (s — 4)/2 of them. 

5.1. Four-Stage, Fourth-Order: RK4()4[4R]. in cases where the number of stages equals the 
number of available storage registers, all possible schemes may be implemented in sR-vdH fashion. For 
the four-stage, fourth-order ER.K, w'e solve lor the minimum error scheme. Setting dA^ fdcn — 0 and 
fU (8) /0ra = 0. results in c 2 « c 3 Rs 4 (5) = 0.0119775. and r T , = 0.613. The relation 


178211381 ’ t3 ~ 222986851 1 ^ . 

between the various other Butcher coefficients may be found in the literature. [9, 21, 35] Stability limits are 
(A, A r ) = (1.42,0.70). Again, a third-order embedded method is impossible with this scheme without FSAL 
constructs. Maximal inner product norm contractivity occurs with Cooper's RK4()4 scheme at c*» & 
and c 3 rs 2^ where r Ti = 1.14373 and A (5) = 0.01755. Gottlieb and Shu[31] use Butcher's[9] case 2 in 
§312, setting b 3 = ^ to. get -1 (5) = 0.01592 and = 0.945. In all RK4()4 cases r^= 0, r £ci = 1, and 
r £ . = 1.393. 


5.2. Five-Stage, Fourth-Order: RK4(3)5[4R+]. An RK4(3)5[4R+] method has 14 DOF, having 
sacrificed only I DOF to low storage. To minimize .1 <5) , Butcher simplifying assumptions ('(2) and 0(1) 
are applied, reducing the constraint, system to 


r'* 0 = 0, A- = 1,2. 3.4, c 5 = 1. r< 5) = e, 6 2 = £f =3 6,c,a, 2 = 0, 

XT/=i ttijCj Cj- /2, / 3,4, — bj[ 1 cy), i = 2,3,4. 
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An exact one-parameter, c 4 , solution has been found where may be made arbitrarily small. For c 0, 
A (5) = ^103/1036800(04 - 1). Unfortunately, both b A and 65 are proportional to (c A - 1)“ J . a so-called 
limiting formula. As c A 1, D = 64 = — [12(c4 — 1)04(604 — 2)] . Setting 5 = —1/40000 and C4 — 199/200, 

RI\4(3)5[4R+]M has A (5) = 0.00003216, A (0) /A (5) ~ 130.3, and D - 6.365. 

To obtain a contractive RK4(3)5[4R+] scheme, we closely follow Kraaijevanger[52] with the excep- 
tion of not enforcing (a 5i - r^a^a^i). Note that we solved 15 equations (8 order conditions and 7 of 
his 8 contractivity conditions) in 15 unknowns whereas Kraaijevanger performed an optimization problem. 
Kraaijevanger's RK4()5 method has = 0.006439, (A,A t ,-) — (1.64.1.34), and (j*£ 2 , r ^oo) 

— (2.191, 1.861, 1.835, 1.508). The RK4(3)5[4R+]N method has A (5) = 0.005635, (A, X v ) = (1.63, 1.40), and 
rj: ) - (1.733, 1.095). We mention that a good embedded method may be added to Kraaijevanger s 
RK4()5 scheme by solving the four third-order embedded order conditions, linear in the 6,'s, by setting 
65 = 113/599. Coefficients and properties of t he two RK4(3)5[4R.+] methods are listed in Tables 3 and 5. 
Stability plots are given in Figure 5.1. 


5.3. Six-Stage, Fifth-Order: RK5(4)G[4R+]. By increasing the stage count to six, a 5(4) pair may 
be considered with Butcher simplifying assumption (7(3). The general R.K5(4)6[4R+] method has 18 DOF 
in the main scheme and 6 DOF in the embedded method. Of the nine main schemes found, the best scheme 
has A (8) = 0.001961 and (A, A v ) = (0.28.0.99). A more ambitious agenda uses only simplifying assumption 
(7(2) while enforcing a condition on To do this we solve 


(5.2) 


r, (fc) = 0, k — 1,2,'*, 5, E U a ^= c </ 2 ' 

b‘> = EEa Mi- 2 - Ei=3 biCidi'2 = EiJ=3 haijOj* = 0, a 42 = 
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A 5 > - 


= 0, 


(3r 4 -rjr~+10r^) 
: 2 (3— I2C3 + IOC3) ’ 


for the main scheme and 
(5.3) 




k = 1,2, 3, 4, 


0, <Z>1 


(«) 


1/130. 


for the embedded method. Note that the strategies described by Papakostas et al.[60] and Hairet et al.[3o] 
(§11.5) must be modified slightly. Tables 3 and 6 show RK5(4)6[4R+]M having A^4 — 0.0009449 and 
(A, X v ) — (0.31,0.93). A stability diagram for this scheme is given in Figure 5.1. With 4>2(? = 1/1440, 
7^ reaches 2 for the RK5(4)6 method. A FSAL method akin to those of Dormand et al.[19, 20] and 
Papakostas et al.[60] is not considered. 

5.4. Seven-Stage. Fifth-Order - RK5(4)7[4R+]. A seven-stage, 5(4) pair may be approached in 
at least four ways: using C(2). (7(3), (7(2) and D(l), or C(3) and D( 1). To satisfy all fifth-order constraints, 
these require 18, 20. 18, and 21 DOF, respectively. For sixth-order these increase to 30, 28, 24, and 25. In 
addition, use of (7(3) reduces the number of embedded order conditions. The simplest approach is to use ( (3) 
and D( 1). Setting ^2x 10“ r> . a solution was found having A* b) = 0.0003974 and (A, A, ) = (0.30, 0.87). 
With only (7(3), a somewhat better solution has A (R) = 0.0003649 and (A, A„) = (0.32.0.89). 

To decrease A (6) further, only C(2) is assumed. Using a FSAL method allows the main scheme to be 
designed independently of the embedded method. For fifth order in the main method, 


(5.4) 


\ k) = 0. k = 1 . 2, 3, 4, 5, 4 4) = = 0, 


b-> = X] i -3 b,(ti2 — X]i = 3 bjCidi-2 — biUijaj-j — 0 , 23j=i a >.i r j 

and for the fourth-order embedded method, 


c“/2, /= 3, 4, 5, 6, 7, 


(5.5) 


r (fc) _ 


= 0, k = 1 , 2, 3, 4, i >2 = ^ Mi*> — r ; ! 


(4) 


0. 
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The remaining degrees of freedom are chosen so that = 0, r f * 6> « -7.8 x Kr 7 , r.*^ ss -I x 10 -, \ 
and ^ 5) = 1/125. The resulting method, RK5(4)8[4R+]FM, has 7 u = 0, A (fi) = 0.00003256, 
- 4<7) = 0.0002906, -4 <8) = 0.0004815, .4< 9 > = 0.0005800. (A. A,) = (0.99,0.98), and (A. A, ) = (1.27,0.81). 
Details of the method are found in Tables 3 and 6, and the stability characteristics are shown in Figure 5.1. 

For phase-lag methods, which we do not pursue, select (<t> 7 — = 0) by setting = 1/756 and 

*4* = 1/7560 [(A, A,,) = (0.36, 1.16)], (a tl = 0 7 = 0) by setting ^ = 1/720 and = 1/5040 [(A, A, ) = 
(0.88, 0.99)], or (a 6 = a 8 = 0) by placing = 1/720 and ^ = 1/5760 [(A,A„) - (0.53, 1.05)]. Note that 
R.K 5 ( 4 ) 8 [4 R +] F M has « 6 as <j> 7 as 0. 


6. Five-Register Schemes. With five registers, the Butcher array is given by 
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and allows (5s — 10) degrees of freedom while having forfeited (s — 4)(,s — 5)/2. 


6.1. Seven- St age. Fifth-Order: RK5(4)7[5R+]. A seven-stage, five- register, 5(4) pair may be 
approached as if it were a 6(4) pair. Both pairs f?(2), D(l) and C ^(3). Z)(l) enable a sixth-order main 
method that requires 24 and 25 DOF, respectively. We will follow the strategy of Sharp and Smart [72] and 
Bogacki and Shampine[6] by solving for the sixth-order method and then will pollute it ever so slightly. For 
a sixth-order main method with C(3) and D( 1) we enforce 


( 6 . 1 ) 


Jfc) 


. =9. k = 1,2, 3, 4, 5, 6, 6 2 = r' 6) = £* =3 fco,a» = £-=3 V?<b 2 

£;=i «ijr] 1 = cf/q, i = 3. 4,5, 6. q - 2, 3, £;* =1 bmj = 5,(1 - cj). 


= 121 j =3 hCidijOj 2 = 0 , 

j = 2, II, 4, 5, 6. cj = 1 . 


and for the fourth-order embedded method. 


rj k) = 0, k — 1, 2,3, 4, 6 2 = £’ =3 &,a <2 = 0, 4><, 5) = 1/125. 

Interestingly, in spite of the nonlinearity in the 5,’s, b 7 = 1/12. Setting r, (6) = 2x 10 -5 , r 8 6) = ?«,••> = 

H2i.j=3 b » c > a ij a j2 -5x10“*, all 20 t} , j) are nearly equally corrupted. The resulting method, RK5(4)7[5R+]M , 
has ,4< fi > = 0.000008959, .4< 7 > = 0.0005771, 3< 8 > = 0.0008997, ,l (y) = 0.001007, A (7 >/A<*> = 64.42, (A,A„) = 
(0.92,0.99), and (A, A,,) = (1.05, 1.19). Fables 3 and 6 and Figure 6.1 display this scheme. 

7. Discussion. In the pursuit of reduced-storage integrators for application to the DNS of compressible 
flow fields, we present 16 different ERK schemes. Schemes vary from third to fifth order in accuracy and 
use from two to five registers of memory per equation per grid point, not including memory used for error 
monitoring/controlling. Schemes have been optimized for accuracy and stability efficiency, linear stability, 
nonlinear stability, dispersion/ dissipation error, error control reliability, and step control stability, all under 
the constraint of reduced memory usage. All presented schemes have been tested by using DETEST, [23] by 
simulating the one-dimensional inviscid wave equation, and by computing standard quantifiable properties of 
tlx 1 Butcher coefficients, as well as using two of the methods in large scale DNS runs. For comparison purposes 
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we have chosen to contrast our third-order schemes to that, of Sharp and Smart[73][SS-R.K3(2)4] . fourth-order 
schemes to that of Prince, [21][P-RK4(3)5], and fifth-order methods to those of Bogacki and Shampine,[6] 
[BS-RK5(4)7], Sharp and Smart, [72] [SS-RK5(4)7]. Dormand and Prince, [19] [DOPRI5-RK5(4)7FM], and 
Papakostas and Papageorgiou,[60] [PP-RK5(4)7F], These reference methods have been chosen because they 
appear to be the best available full-storage methods within their respective classes. The memory requirement 
of these full-storage methods is not less than the stage number for non-FSAL methods or the effective number 
of stages for FSAL methods. 

All schemes presented in this paper have been designed, at a minimum, to avoid any obvious problems. 
As is usual in the design of ER.K methods, great emphasis is placed on reducing A<« +1) to as low as possible. 
DETEST results are well correlated with this measure. DETEST runs involve 25 separate integrations ( Al- 
E5) in 5 general catagories (A-E). Error is computed by t aking the geometric mean of the worst performances 
in each of the 5 catagories by using the Pl-controller. A^ +2) may sometimes be seen to affect, scheme 
performance at lax tolerances. Embedded “quality” parameters B {r+2) , C’ (, ’+-'\ and E (r+ - ] of the low- 
storage schemes are generally quite reasonable, and embedded linear stability domains are commensurate 
with their main methods. The largest Butcher coefficient,, D, never exceeds 7 in any low-storage method 
and for most schemes is near unity. In addition, none of the low-storage methods have defective embedded 
methods. 

Reduced-storage, third-order schemes appear to forfeit little relative to corresponding full-storage schemes. 
At 3 stages, linear stability is identical among all schemes. Accuracy- based efficiency may be brought to 
99% of the maximum achievable with RK3(2)3[2R+]M. Nonlinear stability may be made equal to 84% of 
Fehlberg’s three-stage, third-order method with RK3(2)3[2R+]N while simultaneously requiring 9% less 
work for similar error tolerances. High quality embedded methods are easily added to these schemes. 
Adding a fourth stage to a 3(2) pair appears to lead to 6% higher ?/ (acc) with RK3(2)4[2R+]C relative 
to RK 3 ( 2 ) 3 [2 R +] M . Inviscid stability efficiency also jumps from A/s = 0.290 to A/s = 0.355. If accuracy or 
inviscid stability efficiency is a priority, this scheme is the best, third-order method presented and behaves 
similarly to the 3(2) pair of Sharp and Smart [SS-RK3(2)4]. Efficiencies of these last two methods may be 
seen in Figure 7.1. a comparison of third- and fourth-order methods using DETEST, as well as in Table 4. 
Viscous stability efficiency and contract ivit y, however, favor the three-stage 3(2) pairs. A,,/* = 0.210 versus 
A„/s = 0.175. RK3(2)3[2R+]N has rjr^/s = 0.279, compared to'r^/s = 0.252 for RK3(2)4[2R+]CN, while 
also being 16% more accuracy efficient. Where contract ivity is the primary concern, RK3(2)4[3R+]N nearly 
doubles the normalized contractivity radius of Fehlberg’s RK3(2)3 method (r^/s = 0.333), while still only 
using 3 registers. The price of achieving r^/s = 0.500 is relatively poor »/< acc) , 77% of SS-RK3(2)4. 

A quick survey of existing third-order methods includes several reduced storage met hods by Carpenter 
and Kennedy, [11, 12] Williamson[87], and Wray [90]. Neither the original Williamson nor Wray schemes has 
an embedded method; they have accuracy efficiencies within 0.1% of each other. Of the two methods given 
by Carpenter and Kennedy, both Williamson-type schemes, one is clearly the most accurate third-order 
scheme given in Table 4 but has no error control capabilites, an easily rectifiable matter, while the other 
sacrifices efficiency to achieve an embedded method with no storage penalty. Bogacki and Shampine[5] have 
clearly improved upon Fehlberg’s two 3(2), or 2(3), pairs but the method of Sharp and Smart appears to be 
the best full-storage 3(2) pair. 

Comparing RK4(3)5[2R+]C with the third-order schemes, the fourth-order method is generally not 
only more stability efficient, but a DETEST comparison of all 2R+ methods, given in Figure 7.2, shows 
that it. can achieve moderate error tolerances at, a small fraction of the work needed by the lower order 
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methods. RK4(3)5[2R+]C seems the more prudent choice over any 3(2) pair for all tolerances below fa 10“*. 
Oontractivity aside. RK4(3)5[2R+]C is quite a bargain. 

Optimizing within fourth-order methods may take many directions, with RK4(3)5[2R+]C serving as a 
goo<l reference. Figures 7.1 and 7.3 show DETEST results on the relative efficiencies of all fourth-order 
si hemes and of all three register methods. Table 5 show's that adding a third register, in principle, allows for 
a 6% increase in efficiency with RK4(3)5[3R+]C. Using RK4(3)5[2R+]C or RK4(3)5[3R+]C enables A/s = 
0.334 and A ,,/s ss 0.238. Where accuracy but not stability efficiencies are most important. RK4(3)5[3R+]M 
and RK4(3)5[4R+]M are 22% and 1 1 (5% more efficient according to Table 5. It may be seen in Figure 7. 1 that 
those numbers are not achieved until quite tight tolerances are reached. DETEST results of 4R+ methods, 
Figure 7.4. show that RIv4(3)5[4R+]M is as efficient as RK5(4)6[4R+]M, whose is 02%, to tolerances 
of fa 10“ 8 ! RK4(3)5[4R+]M is acting like a fifth-order method having an »/ acc > of 58% as determined by 
comparing .4 (b >. Both of these 4(3) “M" methods compare favorably with the best contemporary full-storage 
4(3) pair of Prince. [21] Maximum norm oontractivity of fourth-order methods, on a per stage basis, offers 
slightly less possibility than third-order methods. Kraaijevanger’s RK4()5 method is the most contractive 
RK4()5 with r r J s = 0.302, a bit less than Fehlberg’s v r J s = 0.333. This reduction is particularly 
noticeable when additional requirements like low-storage are imposed. With four registers, at least rj /s 
= 0.210 is possible, but this result is likely reduced to r^/s = 0.095 at three registers. These results, along 
with the fact that contractive ERKs do not exist at fifth order, suggest that there is a trade-off between 
oontractivity and order of accuracy. This trade-off may not be so unfortunate because the linear positivity 
radius, r Cx , remains substantial for many high-order methods and it is likely that, the perceived need for 
large revalues is partially attributable to poor temporal error control. Gottlieb and Shu[31] compare two 
second-order methods and find that the noncontractive method, although it has 43.85 times the principal 
error norm of the contractive method performs less well. We inspect existing 4(3) pairs and avoid the 
methods of Fehlberg[26] and Merson[35] because they have defective embedded schemes when used in local 
extrapolation mode. Neither Zonneveld’s method[35] nor Nprsett.’s method[22] are particularly efficient even 
with full storage. The former method may also have an unreliable error estimate on inviseid problems at 
lax tolerances. Even though Stanescu and Habashi[77] offer a 2N method, it lacks both error control and 
efficiency. In the event that overwriting of the ( -vector is not possible, the RK4()5[2N]C method of Carpenter 
and Kennedy, [11] fitted with an embedded method, would be preferable to RK4(3)5[2R+]C because the 2N 
method is 4% more efficient. Compared to Prince's RK4(3)5 method, RK4(3)5[3R+]M is largely the same 
yet. uses only three registers, while RK4(3)5[4R+]M is substantially more efficient . 

The burden of low' storage becomes apparent relative to corresponding contemporary pairs at fifth order 
because of the large number of forsaken degrees of freedom as well as the large amount of research that has 
gone into opt imizing existing 5(4) pairs. This burden may easily be seen in Figure 7.5, a DETEST comparison 
of fifth-order methods. Optimization of lower order methods would seem to have taken a back seat to those 
fifth order and higher for reasons of efficiency. In order to achieve fifth order in 2 registers and 9 stages, 
28 DOF are sacrificed! Not surprisingly, »/ acc ) of 41-45% relative to BS-RK5(4)7 is seen in Table 6. This 
relative inefficiency makes the RK5(4)9[2R+] methods clearly more efficient than the RK4(3)5[2R-t-]C only 
at tolerances of % 10“ 5 to 10“ 6 , and DETEST shows both R.K4(3)5[3R.+]M and RK4(3)5[4R.+]M to always 
be more efficient . To their defense, the RK5(4)9[2R+] methods have been derived with no residual DOF for 
optimizat ion purposes, used no simplifying assumptions, and by virtue of the low-storage strategy, the order 
conditions became horribly nonlinear in the 6, ’s. The brighter side of the relatively high stage number is that 
stability efficiency can be quite high for 5(4) pairs. We hasten to add that if stability efficiency is desired 
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then RR4(3)5[2R+]C should be accurate enough while allowing for much larger time steps. Accepting a 
third register in a fifth-order method allows for accuracy efficiencies to move from 41-45% of BS-RK5(4)7 to 
48-56%, while stability efficiencies stay the same or decline. For acoustic applications, RK5(4)8[3R+]P{8,7} 
ofTers high dispersion and dissipation accuracy on the linear problem while sacrificing nothing on the nonlinear 
problem. When comparing 3R+ schemes, for “M” and "C” methods, fifth-order methods appear to be more 
efficient than fourth-order methods for tolerances less than tn 10 3 to 10 4 . Comparing RK5(4)8[3RT]C 
to the RK5(4)7FC and RK5(4)7FS methods of Dormand and Prince, [19, 20] Table 6 indicates that the 
low-storage method is nearly as accuracy efficient and viscous stability efficient while being more stability- 
efficient on inviscid problems. In this case, the penalty of low-storage is relatively small. One of t he surprises 
in designing low-storage methods was finding b 3 = b 4 = (1 in the RK5(4)8[3R+] methods as well as the 
RI\5(4)7[3R+] method. There are also many other cases of unexpected linear dependencies. We suspect 
that there is an interesting reason behind the order conditions when certain o,- ; — bj, but a theory eludes us. 


Adding a fourth register to a fifth-order method allows for efficiencies that, approach more traditional 
schemes. For RK5(4)6 schemes, RI\5(4)6[4R+]M is arguably better than both of Fehlberg’s methods[26] and 
that of Dormand and Prince[19] in spite of the loss of three DOF to low storage. The most, accurate RK5(4)6 
published seems to be that of Papakostas and Papageorgiou with 4 (6) = 0.0008694, 1.4% better t} (acc) than 
RK5(4)6[4R+]M (.4 (6) = 0.0009449). Sharp[70] offers two RK5(4)6M methods, with A (l,) = 0.0009399 
and ,4 (6) = 0.0009775. He also states that the global minima for RK5(4)6 schemes is .4 (l,) = 0.00087 . 
consistent with what. Papakostas and Papageorgiou have presented. A FSAL method based on RK5(4)6[4R+] 
type schemes has not been pursued. Moving to seven-stage methods, RK5(4)8[4R.+]FM is our only FSAL 
method. With A (6) = 0.00003256, Table 6 suggests that it. is 30%. more accuracy efficient than the DOPRI5. 


Efficiencies based on .4 (7) , ,4 (8) . and ,4 <9) are even more encouraging. The schemes are would be expected 
to perform similarly to compared to Sharp and Smart [SS-RK5(4)7]. Papakostas and Papageorgiou recently 


designed an extremely accurate 5(4) pair [PP-RK5(4)7F] with 6 effective stages. As with the DOPRI5. 
the disadvantage of this approach relative to fully seven-stage methods is the relatively high values of A" 1 
and D , and relatively poor linear stability. On paper, the best 5(4) pair appears to be the Bogacki and 
Shampine [BS-RK5(4)7]. DETEST results show that RK5(4)8[4R+]FM performs as well as or better than 
SS-RK, 5(4)7, PP-RK5(4)7F, DOPRI5, or BS-RK5(4)7 while saving two to three registers of memory. 'Ihese 


results are slightly controller dependent. The threshold for switching from fourth- to fifth-order 4R+ “M 
methods (RK4(3)5[4R+]M and RK5(4)8[4R+]FM) appears to be « 10“ 3 . 


The five-register 5(4) pair RK5(4)7[5R+] is considered to address any »/ (acc) or ,4 (h) shortfall of the 
2R, 3R, and 4R 5(4) pairs relative to existing methods. Designing 5(4) methods based on a sixth-order 
main scheme has been done, first by Sharp and Smart [SS-RK5(4)7] and later by Bogacki and Shampine 
[BS-RK5(4)7], as well as a q(q - 2)-pair by Tsitouras and Papakostas[82] [TP-RI\6(4)7]. For the 5(4) pairs 
4 1 . 1 ’ ) may be set rather arbitrarily, and for these methods .-1 '■ lj ’ is given by 0.9, 7.1. 2.2, and 0.0(x 10 ), 

respectively. What may be a better measure of the accuracy of these methods is -4 (,) . In the same order, 
A< 7 ) for these methods is 5.8, 1.8, 2.1, and 2.1 (x 10 -4 ). Our DETEST results show RK5(4)7[5R+] per- 
forming better than DOPR15, the same as SS-RK5(4)7, and worse than BS-RK5(4)7, PP-R.K5(4)7F. and 
R K 5 ( 4 ) 8 [4 R+] F M . 


It is important to consider the benefits of additional registers so that these benefits may be weighed 
against the cost of the additional memory usage. At. fourth order, switching from RK4(3)5[2R+]C to 
RK4(3)5[3R+]C nets a 6%> efficiency gain. For “M” schemes, RK4(3)5[4R+]M is 126%- more efficient 
than RI\4(3)5[3R+]M in Table 6. Maximum norm contractivity radius increases 130% by going from 
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Rk4(3)5[3R+]N to RK4(3)5[4R+]N, with an attendant 4% loss in accuracy efficiency. With fifth-order 
schemes, moving from RK5(4)9[2R+]( ' to RK5(4)8[3R+]C yields a 12% efficiency gain. Adding registers 
to RK5(4)9[2R.+]M gives a 25% gain with RK5(4)8[3R+]M, 110%. with RK5(4)8[4R+]FM, and 160%. with 
RK5(4)7[5R+]M. 

Below fifth order there does not appear to be a compelling reason to use full-storage methods. At fifth 
order, users must establish the cost of memory relative to CPU time to establish the optimal methods. On 
parallel machines, low-storage methods may enjoy some advantage because of less required communication. 
When sufficient memory is available and fifth-order accuracy is required, RK5(4)8[4R+]FM is essentially 
as good as BS-RK5(4)7, SS-RK5(4)7, PP-RK5(4)7F, and RK5(4)7[5R+]M. Low-storage methods will also 
be relatively more valuable when the number of equations becomes large (i.e. many species). The value 
increases because the storage required of the integrator is directly proportional to the number of integration 
variables yet storage for items like grid metrics is not. 

Stability plots show that step-control stability is enhanced by switching from an I-controller to a PI- 
controller in all of the methods presented as well as the reference methods. Whereas with the I-controller 
schemes are predominantly SC-unstable on their linear stability boundaries, they are predominantly SC- 
stable with the P I-controller. When methods are SC-unstable with the Pl-controiler, it is often at either 
the real axis (viscous) or at the imaginary axis (inviscid), or both. Some room for optimization for each 
of the methods is possible via a and j3. Doing this optimization requires some caution because it is not 
sufficient in the design of a good controller for each of the eigensolutions t.o be damped. The time constants 
associated with these eigensolutions must not be too large or too small. We do not pursue this optimization. 
Possibly a PID-controller could find use in certain DNS runs. Coping with SC-instability is probably best 
accomplished by reducing step sizes. In cases where .4 (<?+2) / J 4 ( ?+ 1 ) > 1, a Pi-controller was found to make 
error control more reliable. Surprisingly, RK4(3)5[4R+]M with ,4 (6) /-4 (5) = 130 was reliable on DETEST 
with both I- and Pi-controllers. In most cases DETEST was able to run at more lax tolerances with the 
Pl-eontroller than the 1-cont.roiler. All low-storage schemes were able to run at tolerances as lax as 10“ 1 
to 10—, except RK4(3)5[3R+]M, which would not run above 10~ 25 with the Pi-controller. BS-RK5(4)7 
had the worst behavior in this regard, possibly because R{z) and R(z) are so similar. With the I-controller. 
DOPRI5. BS-RK5(4)7, and RK5(4)7[5R+]M, especially the last two, had difficulty at lax tolerances. 

Linear advection of information along characteristics is often used as a model problem for studying 
the hyperbolic limit of the Navier-Stokes equations. An extremely difficult test case is the advection of 
information over long distances, because it tests both the spatial and temporal resolving capabilities of a 
scheme. We formulate this test problem with the model equation dU/dt + dU/dx = 0, solved on the interval 
-50 < x < 450. The initial and exact solutions are given by the expression U(x,1) = i exp[-(^=i) 2 ]. The 
exact solution is a wave packet of energy, spread over an interval approximately six units wide, moving with 
unit velocity in time. Note that this test case has information content, at. all wavenumbers. The spatial 
discretization of the first-derivative operator is done with a sixth-order compact operator, known to have 
adequate spat ial resolving capability. The boundary conditions are imposed to ensure that no order reduction 
occurs. [13] 


Figure 7.6 shows linear advection results, obtained with four temporal operators at. three spatial res- 
olutions. The logarithm of the global error is plotted as a function of the work. We assume that the 
spatial resolution dictates the desired accuracy level in the calculation, and that spatial and temporal error 
components should be approximately equal. Note that as the time-step is decreased (increasing work), all 
formulations asymptote to a uniform error that corresponds to the spatial operator component. At, coarse 
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error tolerances (six grid points resolving the wave packet), the CFL condition (temporal stability constraint) 
of all schemes produces temporal and spatial error components that are nearly matched. The fifth-order 
schemes have no apparent advantage over the fourth-order formulations. At moderate and fine error toler- 
ances (12 and 24 four points), the fifth-order formulations become more efficient. The larger CFL condition 
of the fourth-order scheme allows a larger time step, but produces inadequate temporal resolutions. 

To choose a scheme for a DNS run, all of this information must be sorted. First must be established 
the relative cost of memory to CPU time in relation to the CPU and memory requirements of the run. 
The next step is to establish whether the simulation will be more stability bound or accuracy bound. 
Stability bound simulations favor “C” or “S” methods and the 4(3) pairs. For accuracy bound problems, 

methods are probably best, and 5(4) pairs for tighter tolerances. For runs where nonlinear stability is 
deemed important, U N” methods should be used. Acoustic or temporally periodic problems might best use 
“P” methods. Ultimately, // acc > and ?/ stab ^ are the most important quantities. DETEST quantifies r/ acc) 
nicely, independent of order-of-accuracy, while X/s and \ v /s quantify ?; (stab) well. An interesting strategy 
for users may be to choose an acceptable number of registers and then switch between methods of the same 
storage requirements. For instance, at three registers one could use RI\4(3)5[3R.+]C when stability dictates 
the time step and then switch to RK5(4)8[3R+]C as accuracy becomes more important. When accuracy 
is paramount R.K5(4)8[3R-T]M could be used. On stability dominated problems, the general shape of the 
stability domain in terms of (A, A t , ) may be loosely inferred from the stability plots in terms of For the 
sixth-order, tridiagonal derivative operator, the axes on the stability plots may be replaced with <*(~)/2 A 
and — ■ ) )?(c)/4 ^ A t , . This guideline can be misleading at the imaginary axis. W hat tolerance should be 
used for a DNS run? Given the second sentence in the introduction to this paper, atleast 10 -3 would seem 
appropriate. This value also depends on the spatial tolerance, as well as the demands of the phenomena we 
are attempting to resolve. Lax spatial tolerances will negate tight temporal error tolerances. 

It is also useful to consider the effects of simplifying assumptions. Experience in the literature[60] 
suggests that the best schemes are found by using the minimum number of simplifying assumptions. Our 
experience with R.K5(4)6[4R+] and RK5(4)7[4R+] shows that as long as the embedded method can be 
designed, using C( 2) will reduce A {e) substantially over (7(3). Assumption D{\) did not appear to have as 
dramatic an effect. Judging from RK5(4)8[4R.+]FM, it is not unreasonable to think that both BS-RK5(4)7 
and RK5(4)7[5R+]M could be improved upon slightly by using only C( 2). R.K5(4)8[3R+] methods are not 
possible using (7(3). RK5(4)s[2R+] methods have been designed in 9 stages with no simplifying assumptions 
but would require 10 with 0(1) and 12 with (7(2). Adding an extra stage to the minimum number necessary 
for a q(p) pair also appears to be beneficial. [72] 

To demonstrate the usefulness of the methods, both RK4{3)5[2R+]C and RK5(4)9[2R+]S have been 
applied to the DNS of a heated, planar, compressible air jet as well as to methane-air, methanol-air, and 
hydrogen-air flames. \\ 7 e remark that these choices were made long before many of the other schemes 
here were created. In the case of the jet, observing sound generation from the flowfield might be useful. 
Detecting this sound is nontrivial numerically and requires selection of a variable that noticably manifests 
acoustic waves traversing the media. Figure 7.7 shows the volumetric acceleration in this jet flow and the 
sound waves coming off the jet column and leaving the vortical structures. 

Considering an infinitesimal, spherical material volume element , d\ , the volumetric acceleration is given 
by (3 /dr)(D-dr/Dt~) where dr is the infinitesimal radius of the sphere. Figures 7.8 and 7.9 show the 
corresponding vortical and temperature fields. 

An important question in each simulation is at what tolerance does the order-reduction from boundary 
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error show itself. Users seeking tight tolerances would be well advised to consult the literature for known 
solutions to this problem. It may be that hybrid step-controllers similar to the PI- and PID-controllers in 
combination with those for q(q — '2) pairs[8‘2] could add reliability. It would also be very useful to establish the 
stability contours that correspond to rjr x , and rc x , because comparing r ^ to the region of |/f(;)| = 1 
shows t hat J c , is terribly conservative. It grossly underestimates stability on hyperbolic problems. Two of 
these contours would require determining the absolute monotonicity of a polynomial, R{;) or A'(r), with a 
complex argument. 

8. Conclusions. The derivation of low-storage, explicit Runge-Kutta (ERK) schemes has been per- 
formed in t he context of integrat ing t he compressible Navier-Stokes equations via direct numerical simulation 
(DNS). Unlike previous derivations of ERK schemes which focus on only a few characteristics, we attempt to 
optimize methods across a broad range of properties, subject, to varying degrees of memory economization. 
With a storage reduction methodology introduced by van der Houwen and Wray, schemes are optimized for 
stability and accuracy efficiency, linear and nonlinear stability, error control reliability, step change stability, 
and dissipation/dispersion accuracy. The methods in this paper may be reasonably expected to span the 
range of needs for compressible DNS when numerical stiffness is not an issue. 

Sixteen ERK pairs are presented using from tw r o to five registers of memory per equation, per grid 
point, and having accuracies from 3(2) to 5(4). All schemes have high-quality error controllers and generally 
exhibit step change stability when used with a Pi-controller. Methods have been tested by means of not 
only DETEST, but also the ID wave equation. Two of the methods have been applied to the DNS of a 
compressible heated jet as well as methane-air and hydrogen-air flames. Derived 3(2) and 4(3) pairs, where 
few degrees of freedom are sacrificed for low storage, are competitive with existing full-storage methods. 
Generally, 4(3) pairs are more accuracy and stability efficient than 3(2) pairs. W T hen stability efficiency is 
paramount, certain 4(3) pairs are best. For accuracy limited problems, 5(4) pairs are more efficient than 
4(3) pairs as tolerances drop below 10 -3 to 10“ 5 . The transition error tolerance for this switching depends 
on how many registers are being considered. Although a substantial efficiency penalty accompanies use 
of 2R and 3R fifth-order methods because of the enormous number of forfeited degrees of freedom, state- 
of-the-art full-storage methods can be nearly matched while still saving two to three registers of memory. 
Ultimately, the data presented here should help users determine which method is most appropriate based 
on the properties most valued and the relative cost of the CPU time to memory usage. Users will need to 
decide which properties are most valued, make a determination of the relative cost of CPU time to memory, 
and then choose the appropriate method. 
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Appendix A. Implementation of the van der Houwen scheme. 


A.l. Two registers. We now consider the details of implementing a five-stage explicit Runge-Kutta 
method with the van der Houwen methodology for the integration of 

(AD ( JL = F{i jr {t)) 

from time step n to time step n + 1 with only two storage registers. It is understood that / be comprises 
H variables . Third and fourth registers may be used to store an error estimator and the starting V- 
vector. Assume register 1 (fit) contains the /'-vector at time i (n) = / (1) , U {n) = F<0. The function 
1 ■ ^ 1 ) = /' = F* 1 is evaluated and the result is placed in register 2 (fi2). We now perform the 

operations (error estimation and retention of are optional) 

Add = fll 

fie rr = (il -6i)(A<)fi2 
fil = fil + « 21 (A/)fi2 

(A2) fi2 = Rl + (/;, - n 21 )(A/)fi2. 

which translate to 

field = /’ (M) 

fie rr = (6j - b X )(A/)F ( "> 

C (2) = U in > +an l (At)F {n) 

A' (2) = F (2) + {b l - a 21 )(A/)F ( "> 

(All) = //(") +b 1 (At)F< n \ 

where the A -vector is an intermediate vector that is used to pass information from one stage to the next. 
Boundary conditions for the /'('^-vector are evaluated at / (i) = / (n) + c t (A/). This constitutes the end of 
stage 1. 1 he function is now evaluated with the contents of Rl and the result is then overwritten onto Rl. 
With this we compute 

H^rr — Rerv + (6 2 - 6 2 )(A/)fil 

R2= R2 + a 32 (At)R\ 

(A4) fil = fi2 + (6 2 -n 32 )(A/)fil, 


fierr = fi e rr+(6 2 - 6 2 )(A/)F (2) 

= (b-j - 6 2 )(A/)F (2) + (b 1 -bi)(At)F {n) 

F<3) = A -(2) + 032(A<)F (2) 

= F (n) + « 32 (A/)F (2) + bi(At)F {n] 

A -(3) _ V (3) + (6a _ a32 )( A< )F< 2 > 

(A5) = F<") + 6 2 (A/)F (2) + MA/)F (,! >. 
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Stage two is complete. Stage 3 begins with the evaluation of the function with the contents of R.2. Overwriting 
the contents of R. 2 , F -' ! i , with the result of the function evaluation, F- :! \ 


Re rr = Re rr + (*3 - h)(A1)R'2 

Rl = HI + a 43 (At)R2 

(AO) R'2 = Rl + (fc 3 - «. 1;t )(A/)/f2. 

giving 

/ferr = fierr + (<> 3 -fr 3 )(At)F (3) 

= (63 - fc 3 )(Af)F( 3 ) + (6 2 - fe,)(At)F( 2 > + (6, - 6 1 )(At)F (n) 

F (4) = A' (3) + a 43 (A*)F (3) 

= U {n) + a 43 (A/)F (3) + 6 2 (Af)F (2) + 

V (4) _ (jW + (b 3 - a 43 )(At)F {3) 

(A 7 ) = F (n) + 6 3 (A/)F (3) + 6 2 (A/)F (2) + MA/)F <n) . 

To begin stage 4 , the function is now evaluated with the contents of Rl and the result is then overwritten 
into Rl. Hence, 

R err — R err 4 - (64 — 64) (A/)/? 1 
R2 = R2 -f a 5 4( A/)/?l 

(A8) /?! = fl 2 + (64 - a 54 )(A/)/?l, 

or 

Ferr=Fe rr +(t 4 - 64 )(A<)F (4) 

= (64 - &4)(A/)F (4 > + (<»3 - 6 3 )(A<)F (3) 

+ (b 2 - 6 2 )(A 0 F (2) + (/p - 6 1 )(A/)F (n) 

F (5) _ X ( 4 ) + 054 (A/)F (4) 

= F (n) + «54(A/)F (4) + 6 3 (A/)F (3) + 6 2 (A<)F (2) + fe 1 (A/)F ( " ) 

A' (5) = F (5) + (64 - «54)(A<)F (4) 

(A 9 ) = F (r,) + 6 4 (A/)F (4) + MA/)F (3) + MA<)F (2) + 6 ,(A 0 F ( " ) . 

Stage four is finished. On the final stage, stage 5 , the evaluation of the function is done with the contents of 
R.2. Overwriting the contents of R 2 with the result of the function evaluation, we finally arrive at 

Rerr — ^err 4 " (^5 ~ bs)(&t) R2 

(AID) Rl = Rl + 6 5 (A*)fl 2 , 

or 

Rerr = Rerr + ( t 5 - 6 8 )( A<)F <5) 

= (65 - 6 5 )(A 0 F< 5 > + (64 - F|)(A 0 F (4) + (63 - 6 3 )(A<)F <3) 

+ (63 - 6 2 )(A 0 F< 2 > + (bi - & 1 )(A<)F<"> 
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{ r(»+D = v< 5 > + b 5 (At)F {5) 

= l' (n) + b r ,(At)F {5) + MA/)F (4! + b 3 (At)F {3) + b-,{At)F i2) + bi(At)F { ” ) , 

(All) 

where + = /l") -f {At). It may be desirable to write F* n+1 ) back into the register that contained fR " 1 at 

(lie beginning of the time step in cases where the scheme has an even number of stages. If a FSAL scheme 
is being used, then j s use d to compute Fl n + F and 


(A 12) 


Ren-= fien- + {0-be)(At)Rl 


/? err = i? err + ( 0-6 c )(A<)F<” +1 > 

= (0 - 6 s)(A<)F ( ” +1) + (65 - & 5 )(At)F< 5 > + (b 4 ~ M(Af)F < 4 > 
(A 13) + (b 3 - fc 3 )(A/)F (3 > + ( 6 2 - 6 2 )(A<)F (2) + (bj - bj )(At)F <n> . 


Note that register one has F (,,+ 1 1 and that if the step is accepted then F^' + 1 1 = F 1 1 * in the new step. To 
control solution error in a vdH scheme, first some appropriate solution error tolerance is chosen, e 10“ 3 
to 10 -5 . Then one may determine the (A/)l n+1) based on either the I- or PI- step controller. If [f( n+1 ) and 
t 'bi + i) are computed to </ = (p + l)-th and p - th order accuracy, respectively, then we may define at 

time //-hi as + = /J err . Then + is a local truncation error estimate for the lower 

order formula. It is also wise to place a limit on how quickly the time step is allowed to increase, factors of 
between 2 and 5 being the maximum. [69] 

A unique problem of the vdH schemes is that if R a ^ is not employed, then when a step size is taken 
that exceeds the error tolerance it is too late to correct matters. In this case, more conservative values of 
the "safety factor" k might be advised. Normally K = 0.9 is chosen, but this might be reduced slightly here. 
Alternatively, the error tolerance, c, could be reduced so that any transgressions of the reduced tolerance 
might not be a transgression of the original tolerance. It should also be remembered that this procedure 
makes no sense if the F- vector is not normalized in some way so that meaningful comparisons may be made 
between, say, the energy equation and the momentum equations. A possible choice would be 


(A14) 


j(»+i 


rr(H + i) 


in cases where |T r (” + i)| greater than, say, 10 8 (depending on machine precision), and where 
replaces <$* n + 1 > j n £q. (2.15) or (2.16). 

A. 2. Three registers. Extending the vdH concept to allow for three available storage registers for 
a five-stage, non-FSAL ERK scheme, our discussion follows directly from the 2R case but is more terse. 
Assume register 1 (R 1) contains ffb») at time t^ n K The function, is evaluated and the result is placed 

in register 5 (F6). We now perform the operations 


Fold = FI 

R ert = {b i -b l )[At)m 
R1 = Ri + « 21 (A/)F3 


(Air,) 


R2 = m + (b x — «2j )( A/)F3, 
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(A16) 


(All) 


(A18) 


(A19) 

A 


(A20) 


(A21) 


(A22) 


(A23) 


(A24) 


R,,„ = R err + (b-, -!>-,)( At) Rl 
R2 = R'2 + « 32 (A/)7fl + («3i - b l )(At)R3 
773 = 772 + (ft 2 - « 32 )( A/)77J + (ft, - « 31 )(A <)7?3, 

Rx- rr = Rerr + (63 — 63 ) ( A/ ) R'2 

773 = 773 + a 43 (A<)/72 + (« 42 - ft 2 )(A/)771 
R\ = 773 + (63 - « 4 3 )(A /)«2 + {!>■> - a. 12 )(A7)/71, 

« elT = Her,- + (6 4 -6 4 )( A/) 

771 = Rl + «54 ( A/ ) 773 + (053 — ft 3 )(A/)/72 

772 = 771 + (64 - a 5 4 )(A<)ft 3 + (63 - a 53 )(A<)772, 

Re rr = ^err + (65 — 65 ) ( Af ) /7 1 

772 = R'2 + b$(At)R\. 


. 3 . Four registers. 

7?old = 771 

Rerr= ( 6 , -ft,)(A<)/74 
Rl — Rl + « 2 i(A/)/M 
772 = 771 + (61 -a 2 ,)(A/)/M 

/Jerr = Re rr + (6 2 - ft 2 )(A<)771 

772 = 772 + a 32 (A/)/?l + («31 - ft, )( A/) 774 

773 = 7?2+(ft 2 -a 32 )(A/)7?l + (ft 1 — a 3 i )(A7)774 

Rerr — Rerr + (63 63 ) (A/) 772 

773 = 773 + a 43 (A7)7?2 + (« 42 - ft 2 )(A7)7?l + (« 4 i - fti)(A/)774 

774 = R3+ (63 - 043 )(A<)fl 2 + ( 6 2 - « 4 2 )(A <)/?1 + (ft, - a. u )(A()R4 

Re rr = Rerr + (64 — 64 ) ( A M l{'2 

Rl = R4 + « 54 (A/ )/(■:? + (nr , 3 - b 3 )(At)R’2 + (o 52 - ft 2 )(A<)771 
771 = 774 + (ft 4 - a 54 )(A<)7?3 + (63 - a 53 )(A<)fl2 + (ft 2 - « 52 )(A7)771 

77 err = 77 rrr + (65 _ 65 ) (A/) 774 
Rl = 771 + ftr,(A/)7?4 
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A. 4. Five registers. 

Km - K 1 

/ierr = (6, -6,)(AO«5 
m = m + a- n (At)RA 
(A 25) 7?2 = Rl + (6j - a 21 ){At)RA 

H err = R„ r + (fe 2 - b-,){At)IU 
R'2 = R2 + a 3 >(At)Rl + (a 3 i - b 1 )(Al)Rb 
(A 20) RA = R2 + (b-, - o- v ,)(Al)Rl + (6| - a 3l )(At)«5 

Rew = Rfw + (b 3 - b 3 )(At)R2 

RA = RA + a 43 (A/)i?2 + (a 42 - i 2 )(At)«l + (a 41 - h)(At)RA 
(A 27) RA = /?:5 + (i > 3 - a 43 )(Af)fl2 + (6 2 -a 42 )(A*)i?l + {h - a 4 , )(At)/?5 

/ferr — ^err + (^4 — b 4 )(At)R2 

RA = RA + o M (At)«3 + (« 53 - 6 3 )(A<)H2 + (a 52 - b,)(At)R\ + (a 5 , - b l )(At)RA 
RA = RA + (b 4 - a r , A ){At)RA+ (b 3 - a 53 )(At)R2 + (6 2 - a 52 )(At)rtl + (6, - a 8 i)(A<)R5 

(A28) 

Revr — Re rr + (^5 — 6s)(A<)i?4 
(A 29) RA = /?!> + 6s(A<)/?4. 
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Appendix B. Explicit Runge-Kutta Order Conditions. 

Equations of conditions[19] for various orders of accuracy are are found in many places, e.g., §3.4[21]. 
Higher order conditions may be derived by using ButcherMath found in Mathematical##, 89] To provide 
completeness in this work, up to sixth order, these conditions given by 


r i (1) = 

r< 3) = 

t{ 4 > = 

4 4> = 

r?> = 

t (5> ~ 

~3 — 

tP = 

r f> = 

r (») - 
'9 ~ 

r <«> - 

* *) — 

& = 

.! 6) = 

r (6) - 

r 8 “ 

r (R ) _ 
Mo — 

r (6) _ 

- 

r 14 ~ 

r (f>) - 

'16 — 

r (G) ~ 
'18 — 

r <«> _ 

' 20 ~ 


3! 

_1_ 

4! 


_L 

4! 


Ei=l b > 1! 

IE Uibtf- 

I Ei,j = l biOijc) - 

— V s b - — 

5 Eij=i biaijCjdikCk - £ 

1 y 1 * A. a .. f 3_j_ 

6 Z-ij = l u i u *j'-j 5* 

Ei,j,it = l bi a ijCjdjkCk ~ 57 
Ei,j,fc,/ = 1 ^i a ij a jk a kl c l — gT 

i E;',j=i - $ 

6 Ei,j = l b i c i a ij c j — g? 

I E/,j,fc = l b i c h l ij a ikCk 
Et,j,fc = l b i C i a ij C j a jk c k ~ -JT 
I Ef ,j,fe,/= 1 b iaija jk Ck(ljlCl - | 
I E;,M-=1 hidijCjOjk-cl - I 

E,i.U=. b iCi a ij a jk a kl e l — gi 

E;,.,\*.t=i biaija jk c k a kt ci - $ 

Ei,).A-,/,»n=l b i a ij a jk a kl a lm c m ~ 


_(-! 


10 

6 ! 

15 


7v; 


(3) 


rO 


(4) 


^ 4) 

r< 5) 

rf> 

-r 

(0) 

r 3 

r‘ fi > 

r< 6) 

>) 

'9 

_(«) 

”11 

_(«) 

M3 

ws* 

rf ? 1 

J«) 


Ej = l I — tT 
Ei,j = l b i a ij c j — y 
Eij = l biCjdijCj — Jr 
Efj.fc=i b i<*ijajkCk ~ 57 
I Ei,j = l b i c t a u c j ~ f- 
| Eij=l b i c i a i.j c j ~ 5 > 
E(j,*=i hiCidijOjkCk ~ £ 

I E,-j.fe=l hidijdjkcl - £ 

— V' 5 6-r 5 - - 

120 Z-»i = l u * ( * (V 

I Eij,fc=i bjCiajjCjdikCk - § 

1 / 2 111 

2 EiJ.* = l biOijCjOikCk - gr 

X V'-' A. a .. c 1_l 

24 Zji=1 u i u ij<-j fi , 

E s r, . , 10 

ij,kj= 1 biaij a ik<'k<til(‘l “ 

2 E, j,fc=i biOijCjajkCk - g 7 

I Ei=t biCiaijOjkcl - §, 

I E;, , = 1 b > a ij a jk c k - v. 

Ei,j,fc ./=1 b i a i.i c j a jk(lklCl ~ £ 

I Etj,k,i=i bidijdjkakicf - 


V'erner[8G] divides these order conditions into four general categories; quadrature rf k \ k = 1,2,15,4.5.6; 
subquadrature r/, 3) , ly 20 ; extended subquadrature r| 4) , r 2 ^ G 7 , 6 8 10 „ J3 ]4 ir; ]7 18 ; and 

nonlinear x ,. Several higlier-order “tail-tree” conditions of constraint, important in the design of 

linear stability, are given by 


El ^ ^ 0 ij Cl j Ar (Ifc l O { m (l m n C r! _j 

i ,j,A; ,/,m,u= 1 

(8) j i 

7”ll5 — / 0 f ; a jj (I j k^kl^hn^hnn^n o o 

i,j f k y l,m,n,o=\ 

(9) \T^ , i 

^286 — / vj(lj j (Ij k (l k i O-Iji} a ni r( (l no (l 0 p Cp 

i J .k ,1 ,v} ,v ,o,p= 1 
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Table 1: Two-register ERK schemes 



RK3(2)4[2R+]C 

RK4(3)5[2R+]C 

RK5(4)9[2R+]S 

RK5(4)9[2R+]C 

RK5(4)9[2R+]M 


, 1 184746128281*4 

, 970286171893 

i 1107026461565 

i 2756167973529 

, 5573095071601 

tl Jl 

• 3654754301 1857 

I" 4311952581923 

5417078080134 

t 16886029417639 

11304125995793 

a si 

I 394 32254 4 3063 

, 6584761158862 

. 38141181049399 

( 11436141375279 

• 315581365608 

' 70761 55732230 

' 12103376702013 

1 41724347789694 

’ 13592993952163 

1 4729744040249 


346793006927 

1 2251764453980 

1 493273079041 

i 88551658327 

1 8734064225157 

4029903576067 

1 15575788980749 

1 11940823631197 

1 2352971381260 

‘ 30508564569118 



, 26877169314380 

, 1851571280403 

, 1882111988787 

, 6457785058448 


’ 34165994151039 

6147804934346 

“ r 5590444 19395^ 

f 14962850401353 

<*65 



, 11782306865191 

i 846820081679 

j 5771559441664 



1 62590030070788 

1 4754706910573 

1 18167997215013 

a 76 



i 9452544825720 

, 4475289710031 

j 1906712129266 



1 13648368537481 

1 6420120086209 

1 6681214991155 




1 4435885630781 

i 1 18394774831 1 

, 311585568784 




' 26285702406235 

' 9144450320350 

t 2369973437185 

<*98 



1 2357909744247 

1 3307377157135 

4840285693886 



1 11371140753790 

' 131 11544596386 

7758363361725 


, 1017 32471 14-53 

1 1153189308089 

1 2274579626619 

j 1051460336009 

j 549666665015 

° 1 

T 9774461646756 

' 22510343858157 

1 23610510767302 

T 14326298067773 

1 5899839355879 

b-> 

| 8237718656693 

i 1772645290293 

, 693987741272 

, 930517604889 

548816776320 


‘ 136663019*1492 

1 4653164025191 

• 12^94497460941 

’ 7067438519321 

9402906569133 

b 3 

, 57731312506979 

1672844663538 

347131 529483 

31 1910530565 

1 167 2704 94 6363 

f 194048959415^8 

4480602732383 

1509618590291 1 

1 1 769786407153 

1 1 30 L 54 7 166 1974 


101169746363290 

. 21 14624349019 

, 1144057200723 

410144036239 

. 1025420337373 

04 

37734290219643 

1 3568978502595 

1 32081666971 178 

7045999268647 

> 59^0204766762 



1 5198255086312 

, 1562491064753 

1 16692278975653 

. 1524419752016 

05 


• 14908931495163 

' l7 fWfl 74 6 84 f 56 

83604524739127 

» 6755273^90 1 "79 




. 13113619727965 

. 3777666801280 

10259399787359 

Vt. 



' 44346030145118 

' 13181243438959 

43440802207630 

6- 



1 393957816125 

. 286682614203 

j 4242280279850 




■ 7825732611452 

' 12966190094317 

’ 10722460893763 




1 720647959663 

1 3296161604512 

j 1887552771913 




6565743875477 

* 22629905347183 

• 6099056196803 

69 



. 3559252274877 

. 2993490409874 

453873186647 



> 14424734981077 

t 13266828321767 

1528523568003fl 

in 

, 15763415370699 

, 101 6888040809 

i 266888888871 

1 3189770262221 

j 330911065672 

’ 462*70243929542 

'• 7410784769900 

“r 3040372307578 

’ 35077884776239 

' 9937126492277 

b 1 

1 514528521746 

, 11231460423587 

i 34125631160 

I 780043871774 

872991 9304 18 


' 5659431552419 

» 58533540763752 

• 2973680843661 

“f* 11919681558467 

1 1 147305689291 

i, 

< 27030193851939 

1563879915014 

65381 1289250 

483824475979 

i 2575378033706 


' 9429696342944 

68^36 1 07 17585 

9267220972999 

5387739450692 

" l " 14 439^13202205 

b 4 

69544964768955 

. 606302364029 

■ 323544662297 

. 1306553327038 

. 3046892121673 

30 2 C 2026366 149 

t" 971179775848 

' 2461529853637 

' 9528955984871 

1 1013392356255 



1 1097981568119 
r 3980877426909 

j 1105885670474 

| 6521 106697498 

| 1760184658016 



' 4964345317203 

"T" 22565 5 77506855 

’ 89 29 499316 295 




, 1408484642121 

. 14UU555694605 

■ 102651 49U63 

t>c> 



' 8758221613943 

' 19784728594468 

• 209874 J 126425 

i>7 



j 1454774750537 

1 1 183541508418 

. 1643090076625 



' 11112645198328 

» 13436305181271 

’ 4891294^70654 




i 772137014323 

. 3036254792728 

1 116106750067 




~ r 4386814405182 

• 15493572606329 

“r 3955800826265 

h 



1 277420604269 

1 638483435745 

j 866868642257 

w‘j 



1857595682219 

4187244659458 

* 42331321870877 
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Table 2: Three-register ERK schemes 



RK 4 ( 3 ) 5 [3 R+ ] C 

RK4(3)5[3R+]M 

RK4 (3)5 [3R+] N 

RK5(4)8[3R+]e 

RK5(4)8[3R+]P{8,7} 

RK5(4)8[3R+]M 


1 23655924 73904 

, 17396840518954 

1 4745337637855 

, 141236061735 

. 1298271 1761 51 

■ 967290102210 

a 21 

~t" 6146167614645 

1 49788467287365 

' 22386579870409 

T" 3636543850841 

1 60 748409385661 

1 6263494269639 


, 4276267785271 

. 21253110367599 

. 6808157035527 

, 7367658691349 

, 14078610000243 

■ 8 529598215211 

(l 3 2 

6823155464066 

' 14558944785238 

* 13197844641179 

• 25881828075080 

i 41877490110127 

' 5603806251467 



, 4293647616769 

■ 4367509502613 

1 618526949 1 39H 

1 553998884433 

. 8043261511347 

«43 

+ 8986505720531 

1 14519312872408 

1 10454198590847 

1 13597512850793 

1 1150223130613 

1 8583649637008 



8941886866937 

, 1236962429870 

, 2669739616339 

, 15658478150918 

115941139189 

«54 

1 24358012670437 

7464816931160 

t 3429868089329 

1 165836226451 14 

1 92423611770207 






■ 42158992267337 

1 18843935397718 

. 2151445634296 

a 65 

” 


' 

1 9664249073111 

1 7227975568851 

1 7749920058933 





, 970532350048 

• 6206560082614 

, 15619711431787 

a 76 

' 

" 


* 4459675494195 

f 27846110321329 

* 74684159414562 





1 1415616989537 

1 2841125392315 

■ 12444295717863 

a 87 




~T~ 7108576874996 

1 1 484421 76iG0 1 77 

1 11186327299274 



12587430488023 

, 546509042554 

343061 178215 

2491873887327 

, 475331134681 

031 

10870640012513 

11977319897242 

* 9152262712923 

2523 1 50225462 

1 1M 9757507826 

* 7396070923784 



. 6191878339181 

I 625707605167 

4057757969325 

3833614938189 

8677837986029 

a 42 

• 8494387045469 

1 13848262311063 

‘ 5316659119056 

1 8246604264081 

14183712281236 

10519245648862 


■ 3819021186 

. 19121624 16580 1 

1 582400652113 

, 1415160642415 

, 628609886693 

1 2224500752467 

a 53 

2763618202291 

1 12321025968027 

1 7078426004906 

* 13311741862438 

• 81773991 10319 






93461894168145 

4943723744483 

. 1245361422071 

0 64 

” 

" 


25333855312294 

25560747809^6 

‘ 3717267139065 





■ 7285104933991 

, 1024000837540 

■ 1052079198131 

a 7 ^ 

" 

' 


f 14106269434317 

‘ 1998038636351 

> 3768458824028 





4825949463597 

2492809296391 

5225103653628 





1 68 284005 78907 

9064568868273 

856-1 162722535 



, 19773887*15448 

. 314199625218 

■ 514862045033 

■ 346820227625 

■ 83759458317 

b\ 

6523801458457 

1 17714523675943 

• ('19835 0928319 

1 4637360145389 

1 3124407760749 


r 

. 3032295699695 

■ 6528140725453 

1 6410344372641 

o 

0 

0 

l> 2 

' 12397907741132 

' 14879534818174 

' 17000082738695 




r 

| 6 l 26 18 101 729 

, 4395900531415 

1 292278564125 

0 

0 

0 

0.1 

' 6534652265123 

■ 55649460397719 

• 5593752632744 





1 1155491934595 

■ 6567440254656 

, 5010207514426 

o 

0 

0 

64 

> 2954287928812 

' 1.5757960182571 

» 21876007855139 





1 707644755468 

436008689643 

■ 5597675544274 

1 2561084526938 

, 814249513470 

1 6968891091250 

b 5 

' 5028292464395 

9453681332953 

• 18784428342765 

• 79.59061818733 

• 2521463007009 

' 16855527649349 





, 4857652849 

I 195246859987 

, 763521911849 

£>6 

' 

' 


' 7350455163355 

« 15831935944600 

' 8570687289572 





, 10599430 12790 

i 3570596951509 

. 366610483 461.3 

b 7 


" 


' 2822036905401 

9786921605312 

’ 1 1232032898210 





, 2987336121747 

. 1886336382073 

■ 517396786175 

b h 




• 1 56 4565 6 70394 4 

* 9981671730680 

' 6104475350679 

bi 

1 1296459667021 

. 390601394181 

1 1276689330531 

. L 2692994 56 3 1 6 

. 679447319361 

2632078767757 

' 9516889378644 

* 35 030 5155 99 1 6 

' 10575835502045 

' 16631323494719 

' 6240332772531 

9365286548816 

r 

, 2599004989233 

• 31150720071161 

, 267542835879 

0 

0 

0 

u 2 

■ 11990680747819 

’ 68604711794052 

* 12417C7 155676 





1 1882083615375 

, 416927665232 

1 1564039648689 

, 2153976949307 

, 798472430005 

1 138832778584802 

b: 1 

• 8481715831096 

T 69530442^9?41 

' 9024646069760 

' 22364028786708 

> 13882421602211 

’ 30360 4 0.369 75 7 3 

- 

1 1577862909606 

, 3879867616328 

■ 3243722451631 

, 2303038467735 

■ 972791992243 

, 7424139574315 

64 

' 5567358792761 

• 8869216637007 

' 13364844673806 

' 18680122447354 

• 13597677393897 

• 5603229049946 

- 

. 328334985361 

163749046041 

. 606464709716 

. 7354111305649 

• 2994516937385 

3 29 9 3 229351 5 1 5 

*>5 

• 2316973589007 

2599987820560 

• 2447238536635 

’ 15643939971922 

» 00 9 78 5 3 29 5 G 94 

68634 1 5042289 





, 768474111281 

. 1424705874463 

3927384735361 

bn 

" 

' 


< 1008 1205039574 

• 1921 1220871 144 

7982454543710 





1 3439095334143 

, 11199564663291 

1 922429-3159931 

67 

" 

' 


' 10786306938509 

’ 351 36367926059 

’ 15708162311543 





38087261 10015 

1 307718103703 

, 024338737541 

ba 




23644487528593 

13694 1 4400390 l 

' 7691046757191 


35 




Table 3: Four- and five- register ERK schemes 
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Fig. 3.1. Stability limits of the main and embedded methods for' two register schemes , RI\3(2)4[2R+]C, RK4( 3 )5[2R +]C\ 
RK5(4)9[2R+]C, R K 5(4 ) 9[ 2R+] M, RK5(4 )9[2R+]S. Circles denote regions of r^^and r£ 2 contractivity. Shaded regions 
denote locations along the contour |/?(c)j = 1 where the methods art SC-stable with either an f- or Pi-controller . 





Im (z) Im (z) 



FlG. 4.1. Stability limits of the main and embedded methods for three -register schemes, RI\4(3)5[3R+]C , R K 4(3 )5[3R +]M . 
RK4(3)5[3R+]N , /?/\ 5 ;<*/.?/?+/(?, /?A'5(4 )8[3R+]\L RK5(4)8[3R+]P{8, 7}. Circles denote regions of r j ^ , r> 2 . , and 

corifrarfivify. Shaded regions denote locations along the contour |/?(s)| = 1 u'/iere f/it methods are SC-stabk with either 
an l- or P I- controller. 
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Fig. 5.1. Stability limits of the main and embedded methods for four register schemes „ FU\4( 3 )5[4R+]M, RK4( $ )5[4 R +]N* 
RI\5(4)(i[4R+]M, RKi>U)8[4R+]FM. Circles denote regions of r j , ry 2 t ? , £ oc <mrf r/\, amfrar fn'ify. Shaded regions denote 
locations along the contour |/?(c:)| = 1 where the methods are SC-stable with either an I- or Pi-controller. 
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Fig. fi.l. Stability limits of the main and embedded methods for the five-register scheme. RK5(4 )7[5R+]A1. Circles denote 
regions of and r£ , contractivity. Shaded regions denote locations along the contour [/?( z) | = 1 where the methods are 
SC-stable using either an I- or Pi-controller. 
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FlG. 7.2. DETEST Comparison of tu>o~ register schemes. 
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FlG. 7.3. DETEST comparison of three -register schemes. 
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FlO. 7.4. DETEST comparison of four- and fivt-reyisttr schemes. 
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Fig. 7.5. DETEST comparison of 5(4) pairs: low- storage schemes and the reference methods of Sharp and Smart [SS- 
RK5(4)?]' Papakostas and Papage orgmu [PP-RI<5(4)?F], Bogacki and S h ampin e [BS-RK5(4 ) 7]. and Dormand and Prince 
[DOPRI5-RK 5(4 )7FM]. 
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M(». i.H. Ovtrall error on <d />if + Ot'jOx — 0 as a junction of temporal error tolerance, spatial grid density, and nuthod. 
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Fig. 7.7. Volumetric acceleration, 0[ 2 ] = ^ * a ( s ~ 2 )* f or a heated, planar , compressible air jet. The self-similar ,\k jet - 
0.95 air jet exhausts into a quiescent body of air. Sound waves art seen emanating from the jet column and vortical region. 


47 


Meters 



— 22500 

IBs 17500 

1 12500 

7500 

2500 

-2500 

-7500 

-12500 

-17500 

-22500 


Meters 


Fig. t . 8 . \ orticity. V X u (5 * }. for a heated, planar, compressible air jet. 
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FlG. 7.9. Tempt rat urt (°K) for a heated planar, compressible air jet. The T, e t, = 900°K air jet rapidly generates finer 
scale motion due to strong flow instability as it issues into the quiescent T ^ — 300° K air. A htlvin-Helmholt: instability may 
be seen in the downstream portion of what remains of the jet column. 
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Table 4: Third-order ERK schemes: General Characteristics 




Table 5: Fourth-order ERK schemes: General Characteristics 



q o 

o o 

o o 

i- r~ 

o o 

nO CN 
Cl — 

o 

o 

o 

o 

o o 
o o 

o o 
o o 

o o 
o o 

o o 
o o 

o o 
o o 

o o 
o o 

8 8 

o o 

o o 

o o 

rr rr 

o o 

o — 

o 

o 

o o 

o o 

o o 



^ °. 

K K 

o o 

o o 

o o 

o d 

o o 

— — 

o 

o 

o o 

o d 

d o 

o o 

o o 

o c 

ik <k 

q q 




o o' 

X c 

_r 

o 

o o 

o o' 

o’ o’ 

o o 

o o 

o © 

Cl CN 
k, k 

k' tk 

o o 
o o 

(Si x 

o o 
o o 

Cl Ci 
X CN 

o o 
o o 

X Cl 
t- CO 

o 

c\> 

o 

o 

o o 
o o 

o o 
o o 

o q 

o o 
o o 

o o 

o o 

o o 

_; q 

o o 

O — 1 

o d 

- ^ 

o 

o 

o d 

o d 

o d 

o o 

o o 

o o 

8 8 

. * - 

t * co" 


r- — 



u"^ 

1- 

o o 

o’ I" 

CN nO 

rr no 

d d 

o~ o* 







X 

■U3 

o o 

O CO 

rr t- 





1- X 

X rr 

O CO 

(£1 X 

o « 

C: O 

Cl 

CO 

o o 

q rr 

CO X 



cn q 

k 'k 

— — 


— i o 

— -H 

— • o 

— CN 

o 

- 1 

— -1 

— o 

— ' c 

— o 

— ' o 

- 1 






O <K 

ci d 

TT 


X co 

Cl CO 

Cl — 

PO no 

PO r- 

rr Ci 

k Ik 

CN T 
CN CN 

CN CN 
CN X 

O no 
CO X 

ci r~ 
<N M 

— 

•^3 ^}* 

y rr 

t- 

CN 

0-1 

Cl n-s 
CO Cl 

O PO 
CO CN 

O CN 

so CO 


l- X 


CN <N 

CN — i 

~ - 

CN CN 

^ - 

(N CN 


CN 

~ — 


Ci — 



_ 




I 


iT- 


IT- tc 

IT * 

ifi ■£ 

N c 

*n 

IT If 

iTi C 

*r. « 


k k 

k k 

k k 


k k 

k k 

k k 

k k 

k k 

k k 

k k 



* 4 -* 

ua y 

X X 

X X 

X X 

X X 

X X 

X X 

X X 

X X 

X X 

X X 

X X 

X X 

X X 

X X 

h H 






. 




„ , — , 

. — . . 

— . . — . 

■ — • . — 

. — • — 

- — • • — 


CO CO 

re co 

X X 

X X 

no x 




Oj CO 

CO CN 

1 1 

i f 

i i 

1 1 

J J 

1 1 

i i 

1 1 


O <D 


X 1- 

re- 

CO TT 

O Cl 

no" CO 

lO so 


Ci rr 

V k 

r- ci 
rr o 

CO rr 

CO © 

Cl X 

CN PO 

rr o 

no r- 
CN c ;• 

00 Cl 
Cj to 

CN O 
CN 

ce *r 

ce C: 

co 

CO <?! 

Cl 0C 

kD 1- 

co o 

l- r: 

PO Ci 

PO X 

e c 

n : - r \ i 

C-‘- 1— 



^ f-H 

PO Cl 

m cj 

f,*. — 

o 6 d 


rr (N 

t.e 

X X 

— CN 

+ + 

4- 4- 

+ f 

+ + 

■+■ 1 

4- + 

+ + 

4- 4- 

4- 4- 

4- 1 

4- 4- 

4- 4- 

4- 4- 

+ 4- 

cT cm 

+ + 
5 58 

Lk UJ 

o rr 
X X 

no i*0 

rr o 

Co X 
no U0 

X _i 
rr co. 
CN Cl 

O Cj 5 
nQ O 

Ci np 

o o 
o o 

no OC 
rr CO 
co CO 



Cl o 
1 - o 

CN (N 

X PO 

rr PO 

no p‘i 
PO PO 

O 1- 

no t- 
o — 1 
q po 

o r- 

O X 

— o 


o o 

o d 

o d 

O O 

o o 

o o 



o o 

o d 

— o 

CN — ' 

o o 









X 

t- 

f- 

o 

to 

o 

ce 








Ci 


CN 

l n 

o 





Q 

X 

o 

nO 

X 

00 

OC 

— 

O 

CN 

o 

tO 





o 



o 

CO 

o 


- 

“ 

ilO 














CO o 

X o 

X' X 

t- PO 

X o 


4- 4- 

5 58 

X t- 
Oi X 

© co 
X no 

d cv 

1- CO 
Ci X 

to o» 

Ci — 
X Ci 

. . 

■ ■ 

V O 

-4 o 

l~ o 

X X 

O CN 

q ^ 

r- q 


o o 

o o 

d o 

d o 

— ^ 

o o 



— 

o o 

— o 

—1 o 

— — 


\J 'sj 















r i .?i 






f- o 



rr o 

— o 

X o 


1- 


4- + 

O CN 
rf O 

nO O 

l~ CN 
01 Cj 

rr 6 

IaO 8N 

ce- ce* 

O Ci 
Pi Cl 


- - 

P'.| o 
o 00 

<N O 
Ci X 

x *r 

OO o 

X PO 


S_ _ P 

CD C5 

-H — 

— t O 

o o 

— 

-H 

— d 



— 1 o 

o d 


CN ci 


1 


co o: 

X* X 

X X 

ce* c’v 

CO CO 

XX 



CN CN 

| | 

(N CN 

| | 

0; Cj 

| | 

X PO 

1 1 

7 i 


+ + 

«-T* <rr 

— i <n 
x no 

X' CO 

CN (N 

ci co 

rr nO 

o o 

1- X 

co co 

O QD 

^ TT 
O 0C' 

o: o 

tm 

p: ii 

r- a 



o r- 

X o 

rr o 

tf u i 

PO <N 

X CN 
Cl no 
X O 

X X 
n'' X 

PO o 

rr no 
rr CN 

H 


X no 


i- -*r 

TT 

o; o: 

00 no 



, *v 


(ii no 

rr 

—» ^ 



o 

o 

o 

o 

o 

O 

• 

• 

o 

o 

- 

- 

o 


i 

X. ^ 
' CN 
, CN 

~ H C 

CN 
Ci 2 
° o 

0,68 

0.136 

n*> ® 
” CN 

"" O 

£ 3 

° d 

n'l ° 
— O 

- ■ 

• • 

0.63 

0.126 

0.62 

0.124 

no © 

x h 
° o 

x 3 

x q 
° o 

1- z 

° o 

0.81 

0.162 

<< — 

© £ 
io 

. - o 

i- £ 

__ X 

ZL x 

r>-. ^0 

® S5 

5 s 

. o 

<£> |J 



o© ^ 

co 2 

£ 8 
• o 

— CN 
„ CN 

- o 

X £ 

° o 

a ; 

- o 

<-< 

^ c 

" -l o 

° d 

— o 

° o 

^ 0 



° o 

° o 

° o 

° o 

° d 

■ 


X 





CN 


Cv 

« 

o 

no 

X 

CN 






tf 

OC 


to 


Cj 


r- 




« r 

X ' 

X 1 

c " 

°0 4 


X • 

n 1 

00 1 

CO ’ 

o ■ 

r- * 

1- • 

X 1 

LJBmI 


o 

d 


6 

tN 

o 

o 

o 

o 


o 

o 

o 

Hi 


•CO CO 


cn" X 

X CO' 

X X 

XX 

XX 

X X 

CN PO 

PO X 

X X 

X PO 

CN P0 



1 1 

1 \ 

1 1 

1 1 

1 1 

i I 

1 f 

1 1 

1 1 

1 1 

1 1 

1 1 

1 1 


4 - 4 - 

5 5 Hi 

1 - o 

- o’ 

00 o 

I- CN 

Ci CN 

o o 

Tf H 

Cl Cl 
1- 

o p: 

QO TT 

r- tt 

CN Cl 
CN X 

nO O 
rr CO 

>0 t- 
— ' X' 

rr rr 
O rr 

<01 no 
Ci CN 


-*r x 

1- nO 

X no 

o t- 

— i 1- 

Cl ^ 

rr Cl 

C1 rr 

X Cl 


rr Cl 

Ci rr 

X Ci 



1- CN 

rr ci 

^ 00 

up CN 

rr CN 

^L> Oj 

PO -H 

rr <oi 

—i d 

PO •-* 

I- PO 

Cl rr 

— X 



ce> co 

XX 

re re 

x - xT 

i i 

i r 

XX 

X PC' 

esi co 

co ro 

| | 

X PO 

1 1 

CO co 

1 1 

CN PO 

1 1 


+ + 

5 58 

1 1 

— ’ cT 

(N I-. 

1 1 

co X 
no X 

TT O 

1 1 
CO CN 

1 1 

CO CN 
Ci 

1 1 

O to 
ro oo 

^ CO 

•V o 

o’ X 
O pj 

t- to 
O re 

x r- 

C1 CO 

CM no 

rr <£> 

ce. i- 

o’ xT 
no P* 


-r -r 

^ TT 

X X 

CO Ci 


CN CO 

o o 

X no 

CN X 

(N X 

a- co 

-T X 

^ r *" 

rr PO 



to ce« 

ro x 

- ^ 

^ X 

PO -H 

o co 

p: cn 

rr X 

— X 

— . -H 

1- rr 

Cl no 

-r X 

1 

* 

_ CN 
CN 

-i 

. CN 

« 2 

s s 

' , C'l 

« ® 

a ° 
£ « 
T CN 

i£ 

«■ 

i- £ 
CN 

e! 

_ CN 
X 50 

pS o 

P CN 

_ Cl 
X £ 

o § 




^ o 

° O 

■“ o 

° o 

“ o 

° o 

*“ o 

° o 

° o 

o 

° o 

° o 

Hi 

-< "d 

1- ' r 

S3 

1- 5 
x a,-- 

CN J 
^ O 

V. i 

CN ^ 

^ o 

« CN 
. P0 

p*. °5 
oo 

zt 

« X 

si 

CN ^ 

CN 
CO $$ 

PO g 
no 

hi cc 

” CM 

B 

•K 

o 

■“' d 

° o 

^ d 

° o 

^ o 

° o 

”* o 

d 

° o 

^ o 

- -1 o 

d 

1 


0 

‘vj 

y. 











1 

\* 

4- 

+ * 


+ 

+ 

+ 

S' 

S' 

no 


iO 

t»T 

kl" 


Ik 

a; 

cn 

n0 

ft; 

•ro 

to 

a; 

X 

c; 

tO 

a; 

V 

C5 

no 

CN 

CO 

CN 

nO 

■'? 

0‘ 

ce- 

TT 

re 

*r 


o* 

CO 


— 

- 

rr 

^r 

TP 

rr - 


y 












VJ 


tt; 

»--« 

ft; 

a; 

a; 


cq 

rr 
' > 


■*r 


^r 


CD 

0D 


K, 






s 

tt 

5 

S 

a? 

5 























■ 













CN 










<L> 

1/ 





-T- 









d 

<L 



X' 

bC 

"ai 

















Author 

Year 

c 

© 

Du 

c 

u 

ol 

Ck 

c 

o 

i> 

CL 

c 

Si 

Wj 

4i 

Du 

c 

ai 

Cl 

c 

u 

CL 

C co 

rj Ci 
Cl 

m — 

E 

a. « 

S 05 
u 2 

L 

2 «o 
1 2 

Prince 

1979 

JO 

3 S 

b! 2 

-C 

— CO 

x 
t Cl 

Du — 

c 

S X 

s ^ 

O Ci 

N - 

o 

'£ 

„ no 

S 2 


51 


Table 6: Fifth-order schemes: General Characteristic 
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